a 


For Reference a 
NOT TO BE TAKEN FROM THIS ROOM nnaene 4 if tf 


M 4 
Late ingeaaty 
igvars hens aie bsiranguah we vi 
SPR OaR URS AD eietined dik 
Uelihcier paa 


Wawona, 
NG ae 


Paces 


Spa atys gett eas pa kane oug NMsiEgIAG es 


Bectrts pasys 
aiesciirs 


Digitized by the Internet Archive 
In 2024 with funding from 
University of Alberta Library 


https://archive.org/details/Milroy1976 


Tae ey Us Ne Tea hath ene ens Ce AGB BR a A 


RELEASE FORM 


NAME OF AUTHOR Richard Douglas Milroy 

TITLE OF THESIS Beat Frequency Heating of a Plasma Near 
the Electron Plasma Frequency 

DEGREE FOR WHICH THESIS WAS PRESENTED Master of Science 


YEAR THIS DEGREE GRANTED 1976 


Permission is hereby granted to THE UNIVERSITY OF 
ALBERTA LIBRARY to reproduce single copies of this 
thesis and to lend or sell such copies for private, 


scholarly or scientific research purposes only. 


The author reserves other publication rights, and 
neither the thesis nor extensive extracts from it may 
be printed or otherwise reproduced without the 


author's written permission. 


er a a Aap ahires \reecce*t eek 
grVicAe-, een’) wepreew) «4 

op it=9 Sy (Gee seaside ve Sey red-oak 
eer Aneveers Weamed ore 


jreteaiot weds Paladvgicthim & o-taacesgt 
uae fa Mao bist ipempaiper s¢ fh. t4tQ eee. | 


his bey ~ cigs teaw Bie =o heét of tes ee 


2 Ww 


a 


THE UNIVERSITY OF ALBERTA 


BEAT FREQUENCY HEATING OF A PLASMA 


NEAR THE ELECTRON PLASMA FREQUENCY 


by 


RICHARD DOUGLAS MILROY 


A THESIS 
SUBMITTED TO THE FACULTY OF GRADUATE STUDIES AND RESEARCH 
IN PARTIAL FULFILMENT OF THE REQUIREMENTS FOR THE DEGREE OF 


MASTER OF SCIENCE 


DEPARTMENT OF ELECTRICAL ENGINEERING 


EDMONTON, ALBERTA 


SPRING, 1976 


sDarogaye UNESRT FONTINGS ENT GARR, 


ad 


FD eatie acide, abies 


aunwws @ 

ipusheae Sep THOME CrANGARD %) Trsypet FHT OF: 

40-5396 THT aot SewseatiwoEs ape fo tesR.TISzeF. 
Ateatou Yo pRee oe 


QATANTE ERED LEO M4ASt WH Pere" tites 


>. 


THE UNIVERSITY OF ALBERTA 


FACULTY OF GRADUATE STUDIES AND RESEARCH 


The undersigned certify that they have read, and 
recommend to the Faculty of Graduate Studies and Research, 
for acceptance, a thesis entitled "Beat Frequency Heating of 
a Plasma Near the Electron Plasma Frequency," submitted by 
Richard Douglas Milroy in partial fulfilment of the 
requirements for the degree of Master of Science in 


Electrical Engineering. 


naa’ anba \adae 19d Sena elidw> tarrbemetee eff 
phaareitt Sas peknaid? sdeuteks 40 otter acd of teenT 
Pe Fiptes parngedss Tat? tol iicas chases’ » ,o0es 190088 TERS 7 
42 corlgel. © CQ btieeee: Ge80l4 sostross ad? 2098 seeale @ 
® 


IS: “| #6950 selquvod 


- 2 
A 


ABSTRACT 


The heating of a plasma in a solenoid, with a beat 
frequency harmonic which is excited at a frequency near to 
that of a Langmuir mode in a plasma, is examined. m3 Gus 
shown that at high temperatures the heating rate is very 
insensitive to changes in plasma density. The amount of 
energy that can be coupled to a plasma in a solenoid with 
this heating scheme is investigated by using a one 
dimensional computer code which incorporates an exact 
solution of the relevent MHD equations. The absorption of 
energy from a high powered laser is shown to be 


Significantly enhanced with this process. 
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CHAPTER 1 


INTRODUCTION 


The heating of plasmas to thermonuclear temperatures is 
a major problem in the quest for controlled thermonuclear 
fusion. A problem exists with many of the heating schemes 
that have been used in experiments in that their efficiency 
drops off with increasing temperature. Plasmas in a 
solenoid are often heated through inverse bremsstrahlung 
absorption of laser radiation. It is difficult to achieve 
good energy coupling at high temperatures with this 
mechanism, since the absorption length is proportional to 
72 


T,'s In this thesis a new technique of heating a plasma in a 


solenoid is investigated. 


This new technique involves enhanced absorption of 
radiation by a plasma and can be achieved when a difference 
frequency harmonic is generated through the nonlinear mixing 
of two intense laser beams. The power absorption can be 
optimized by forcing a mode which is close to a natural mode 
in the plasma. In this thesis the case of two antiparallel 
lasers with a difference frequency close to that of the 


natural Langmuir mode is investigated. 


Theoretical studies of heating through the nonlinear 
mixing of waves close to the plasma frequency have been done 
by Fuchs et al.!, Cohen?, Schmidt3, and Capjack and James*. 
The results of Capjack and James allow a heating rate to be 


calculated from the laser powers, the plasma temperature, 
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and the plasma density. 


Beat frequency heating could prove useful as a 
technique in heating plasmas to very high temperatures. 
This approach is in contrast to the normal heating approach 
where inverse bremsstrahlung is the only heating mechanisn. 
The absorption length of laser radiation through inverse 
bremsstrahlung is very long at high temperatures, so the 
path length of the laser beam through the plasma must also 
be very long to ensure that a reasonable fraction of the 
laser energy is deposited in the plasma. The nonlinear 
heating process, however, could allow a significant fraction 
of the laser energy to be deposited in a short length of 


plasma. 


Studies to date on the nonlinear heating process assume 
that the frequency of the difference frequency harmonic is 
very close to the frequency of the Langmuir mode. In most 
experimental situations the frequency of the Langmuir mode 
will change significantly during the heating process due to 
a change in electron temperature and density. For this 
reason it impossible for the frequency of the difference 
frequency harmonic to be close to the frequency of the 
Langmuir mode during the entire heating process. on this 
thesis an investigation is made of the sensitivity of the 
nonlinear heating process to changes in density. The 
coupling of energy to a plasma will be shown to be increased 


by increasing the plasma density beyond the value necessary 
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to make the frequency of the Langmuir mode equal to the 


difference frequency of the lasers. 


The amount of energy that can be absorbed by a plasma 
in a solenoid through the nonlinear process is compared to 
the amount of energy that can be absorbed through inverse 
bremsstrahlung. A computer code which Simulates the 
magnetohydrodynamic behavior of the plasma as it expands 
outwards during the heating process is used. The Simulation 
allows an assessment to be made of the plasma conditions 


necessary to optimize the heating process. 
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CHART ERe 


THE HEATING MECHANISM 


2.1 PHYSICAL MECHANISM 


In this section some of the results of a paper by 
Capjack and James* on the use of beat frequency heating at 


the electron plasma frequency are briefly summarized. 


In their study the following assumptions are made: the 
heating due to the nonlinear interaction of two high 
frequency laser beams is considered; the source waves have a 
frequency much greater than the electron plasma frequency 
and so they can be approximated as free space waves in the 
plasma; the electric fields of the source waves are linearly 
polarized with components only in the y-direction; the dc 
magnetic field B is in the z-direction; the plasma has an 


isotropic Maxwellian equilibrium distribution function. 


The analysis begins by finding the plasma distribution 


function from the Vlasov equation. 


ote Bf eee VXB aro tee 
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where e=+ for ions and - for electrons. The solution to 


equation (2.1) is expanded as 
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Here the frequencies and propagation vectors of the two 
source lasers are (w, 7k,) and (woek,)- The frequency and 


propagation vector of the driven beat frequency wave are 


(w3¢k3) where 


Dae OT =e 


k3 = ki - ke 


(2. 3) 


It is assumed that Boe kiv Kor and kk3 are all in the 


z-direction. 


An expression for the zero velocity moment of f3° is 


found to be 


(2.4) 
60 
aes, Ei Cavan, Up ias Cu) 
+ ? (— fe) Ae boa du 
V 9) (w3 = K3 u) 
iS W] WO be Be 
where u is the velocity component along the z-axis and ie 


is the thermal velocity. If this result is substituted into 


Poisson's equation 
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3 fos 3 
V+ E; = 4mie {{f3 (v) dv - Jf3 (v) dv} (2.5) 


and the nonlinear term for ions is neglected, the following 


expression is found for F3 . 
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where 
n. = plasma density in cm-3, 


0 

K=k3A 5 where hy = Debye length, 

A, and 4, are the wavelengths of the incident laser beams in 
microns, 


P, and P, are powers of incident beams in Watts/cmn@, 
ao. = w3/(k3Vg,) 
sale eS 
iks © expl-m, (U4 ky Te 


Fo =— : 
eee ae (w3 - k3 u) eta 


StixS has shown that Ee can be written as 
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Fo YT rep expl-ag, 1+ 2 1 Slay.) (2-8) 
where 
(2.9) 
The function S(z) is just the complex error function. 


Collisional damping effects are accounted for by defining 


the real part of Pe as 


d > 
Ro(Fy ) = Max{|—="—, exp[-a_ ] (2. 10) 
2w 30 |k3| ee 
ie) 
The electron-ion collision frequency, Yas may be 
approximated as (see Clavier®) 
-5 
BO SK 10 n (2.11) 
at hia Siete Jo) a) 
el a" 


e 
where the plasma density n), is given in cm-3 and T, in 


electron volts. 


The denominator of equation (2.6) can be written as 
K2p, where D(k3,w) is as defined in Krall and Trivelpiece’, 


equation (8.3.11), as 
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D(k3,w) is called the dielectric function, and its zeros 
correspond to natural modes of the plasma. If a Maxwellian 
distribution function is substituted into equation (2.12) 
for af and the integration is carried out it can be shown 


that 


?2 2 
KD=K +2+i(a. FF.” + an (2. 13) 


The- average power absorbed by the electrons in the 


plasma from the induced signal is given by 


ita 
W3 = Ra {(J3z) E3} (2.14) 


where J3z is the induced electron current in the z-direction 


(see discusion in Capjack and James® following equation 


(11)). 


By expressing J3z in terms of the mobility tensor as 


given in Stix5, the heating rate may be expressed as 


2 
a 5 
e - 
Ws = b= R, (Mz) IEs| (2.15) 
where w is the electron plasma frequency, %- is the 


a. a. 7 oe Danan 
ee a a 
vdteas Gir 405 nokta aizapeteit | a4? bebtes al | 


setilvried » 21 al edt 25 e0hos [etefee Gt feo: 


j 
‘0 


ve 


PQ) amiteins wid Geen dicdict «1 aal see? paites | 


29062 88 a> +1 tua tity 6B odsaayese! 029 bas y 


. «£ 


ite ° ® yj 
oem * a) ugh f esr “Ye 0i 
_ 
2 ni 60) San] wie we h1980S Jepe4 922648 BAT 
ye eels at fergée S42000: add coed Gee 
F — : 
7 ' 
er.) : — a 
igs (seb)? ai* € 
© af 
a 
ih 44 =~ 440% 7 = : 4 
Sher oer bet ty Saat aes ates “SS DOWERL off ai 5% * 


Ho peioe orlesolint 49eee6 beg aserqas al 
——— 
- ee Pr 


colenoesan @ 


- 


oe Tones tratidos mae 7” mate: wilh, 


_ 


=e tngeagte : 


nits 
Bh. Njeeprare ‘at You. aie ted mtr trite. 


+h ati r re a wir 
_ a a J 


- 
: 7 - 


i a 


het > Paes 
- : 
So * A iostosia = ie 


ae iu 


electron cyclotron frequency and hee is the electron 


mobility tensor. Substitution for ae yields 


ee 
= O, ee k 2 
_ pe o- ei 3 -0, 2 ergs 
W3 Wena Max} ———- 3 vn e O-| |E,| (37> —) (2-16) 
beam 2 30, [ks] cm sec 


This equation is valid in the limit where either collisional 
or collisionless damping is dominant. The expression is 
only an approximation in the region where collisional and 


collisionless damping effects are of the same order. 


PEP HEATING CHARACTERISTICS AND 


THE LANDAU DISPERSION RELATION 


When the central core of a plasma in a solenoid is 
heated, the plasma expands outwards leaving a reduced plasma 
density in the core. The rate at which the plasma expands 
outwards during the heating is dependent on the heating rate 
and the initial magnetic field in the solenoid. Situations 
Can exist where the density of the central core is depleted 
to less than ten percent of its initial value. Thus, since 
the density and temperature of the plasma can vary greatly 
during the heating process, the frequency of the Langmuir 
mode can also vary greatly. Before doing a Simulation of 
such an experiment, a knowledge of the dependence of the 
nonlinear heating mechanism on the plasma density and 


temperature is desirable. This allows the initial plasma 
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parameters to be chosen so as to optimize the plasma heating 


rate. 


The heating rates can be calculated by numerically 
evaluating equations (2.6) and (2.14). The specific case 
that will be considered is one where two CO,, laser beams 
with wavelengths of 9.6 and 10.6 microns are mixed in a 


hydrogen plasma. The heating rates calculated below are 


given in terms of the dimensionless quantity 


peat Zed 
W3 = G) ALE (2.17) 
wherey,; = 1 sec. watts? /(ergs cm*). W, is proportional to 


the power absorbed/particle in ergs/sec if P,; and P. are in 
watts/cm2. The electric field is normalized so that the 


quantity, €,, is independent of the incident laser powers. 


E3 (2.18) 
VP; Po 


E3°= 


The graph in Fig. 1 gives the values of W 3 and £3 
versus n for various values of temperature, T(eV). The 
graph in Fig. 2 gives the values of W3; and E3 versus T(eV) 


for various values of density, no: 


The sharp minimum at about 50 eV in the second set of 


graphs corresponds to a point of minimum damping of the 


@riven wave. Here collisional effects are becoming small 


and collisionless damping is just taking over. 
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(solid line) and difference frequency 
electric field (dotted line) versus 
temperature. The plasma density is 
indicated on the curves. The case where 
n, =1.08x1017 cm—17 corresponds to 
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From equation (2.16) it may be seen that the rate of 
energy coupling to the plasma is maximized when 1E31/2 is 
maximized through the selection of an operating point. With 


the use of equation (2.13), equation (2.6) can be written as 


VP PS 


1 3.11x107! kgayag(14i a FO) 


ic ty 
K D(K,ws) 

Thus {E31!2 is maximized if K2D is minimized. Minimum values 
of K2D are expected when the driving frequency, W; , is close 
to a zero of D, or in other words, when the the plasma is 
forced at a frequency close to a natural mode in the plasma. 
The position of the zeros of D(K,w) is a function of K. The 
parameter, K, is a function of density and temperature and 


can be written as 


cael eee ee 
Reeekat ht) [es rer nag 
D 4re 0) y2 (wo/ks) 


This suggests that the response of the heating process to 
varying parameters such as density and temperature can he 


qualitatively determined from a plot of the zeros of D. 


An iterative procedure used by Capjack and James* was 
used to find the first three zeros of this expression over a 
tange of K values. The first zero of this expression yields 
the complex frequency corresponding to the Langmuir mode. 


The other zeroes correspond to the complex frequencies of 


plasma oscillations which are not normally considered as 


natural modes of oscillation since they damp out in 
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approximately one period of oscillation. 


Figs. 3(a) and 3(b) are plots of complex frequencies 
corresponding to the first three roots of D(K,w), normalized 
to w3 4 aS a function of density, for temperatures of 100 eV 
and 500 eV. Figs. 4(a) to 4(c) are similar plots, but are 
plotted as a function of temperature for plasma densities of 
4x102© cm-3, 1.08x10!17 cm-3, and 4x10!7 cm-3. The "x! at 
1.0 on the real axis of these plots indicates the position 


of the pole corresponding to the driving frequency, wu; .- 


Fig. 1 shows a sharp resonance in —; and W3; for a 
density of about 9x10!© cm-3 when the electron temperature 
is 100 .eV. Fig. 3 (a) shows that this corresponds to the 
point where wz; is closest to the natural frequency of the 
Langmuir mode. Fig. 1 also shows that no such resonance 
exists for temperatures higher than 500 eV. From Fig. 3(b) 
it can be seen that the reason for this is that the Langmuir 
mode is heavily damped at densities that make w; equal to 
the real part of the Langmuir mode's frequency, and thus its 


pole never gets close to w3- 
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FIG. 3 (a) Frequency and damping rates 
of natural modes near the electron 
plasma frequency 4S a function of 
density. The temperature is 100 eV. 
The complex frequencies of the natural 
modes in the plasma are given by 
w = w3 (é-i iy). Value of parameter 6 is 
given in brackets where 6 is defined so 
that the density, 0, =§x1.07x10!7. 
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FIG. 3(b) Frequency and damping rates 
of natural modes near the electron 
plasma frequency as a function of 
density. The temperature is 500 eV. 
The complex frequencies of the natural 
modes in the plasma are given by 
w = w3( iY); Value of parameter 6 is 
given in brackets where 6 is defined so 
that the density, esate ca NUE ABER 
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FIG. 4 (a) Frequency and damping rates 
of natural modes near the electron 
plasma frequency, aS a function of 
temperature. The density is 4.0x101!6 
cm-3. The complex frequencies of the 
natural modes in the plasma are given by 
w= w3(é-iy) - Value of temperature in 
(ev) is given in brackets. 


17 


- 
_ 2 ; 7 eet 7 — A 

Per aes) 924 iirpeTd fs 

7 i aS SoG, leapt 


3-00 4e50 6-00 JEoD 


-50 


Mode 2. 


Langmuir Mode 
(3000) & “4 
(2000) Pr 


(1000) g 


FIG. 4(b) Frequency and damping rates 
of natural modes near the electron 
plasma frequency, as a function of 
temperature. The density is 1.08x10!? 
cm—-3. The complex frequencies of the 
natural modes in the plasma are given by 
w= w3(é-iy). Value of temperature in 
(eV) is given in brackets. 
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FIG. 4(c) Frequency and damping rates 
of natural modes near the electron 
plasma frequency, aS a function of 
temperature. The density is 4.0x10!? 
cm-3. The complex frequencies of the 
natural modes in the plasma are given by 
w= w3(E-1iy)- Value of temperature in 
(eV) is given in brackets. 
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In Fig. 2 it can be seen that there is a resonance in 
E3; for a temperature between 100 and 200 eV and a density of 
4x102© cm-~3. Fig. 4(a) shows that the distance between the 
pole at w3 and the pole corresponding to the Langmuir mode 
is minimized in this temperature range. The general shape 
of the graph of £3; versus temperature for a density of 
1.08x10!7 cm~3 can be explained with the help of Fig. 4(b). 
The increase in —; as temperature increases to 6 eV is due 
to a decrease in tollisional damping. The decrease in E3; as 
temperature increases from that point can be attributed to 
the fact that the pole corresponding to the Langmuir mode is 
moving further from w3 with increasing temperature. Hig. 
shows that at a density of 4x10!7 cm-3, £3; is almost 
temperature independent from 1 to 1000 eV and then starts to 
drop off. Fig. 4(c) shows this to be consistent with the 
motion of the Langmuir pole away fromw3 - Fige 4(c) shows, 
however, that at low temperatures poles corresponding to 
other natural modes of the plasma are close to w;- In fact, 
at a temperature of about 50 eV, the distance between w3 and 
the pole corresponding to the second mode is minimized. As 
the temperature increases this pole moveS awav from w3 
rapidly. Thus, “some” "sort Ylof” .resondnce Sinm €,eetoreead 
temperature close to 100 eV would be expected. THISSTS?) “Now 
observed. Because of this discrepancy the importance of 
these higher order modes in the heating process was 
investigated further. Contour maps of the function 


{K2D(K,w) | for the quadrant ke (w)>O and Im(w) <0 are shown in 
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Figs. 5(a) to 5(e). These contours are drawn in such a way 
that |{K*D(K,w) |) changes by equal amounts between contour 
lines. A set of computer generated three dimensional plots 


of 1/(1K2D(K,w)|) is shown in Figs. 6 (a) to 6(e). 


Figs. 5(a) - 5(c) and 6(a) - 6{c) show that for large 
K, higher order poles are almost as wide as the pole 
corresponding to the Langmuir mode, but they are heavily 
damped. It is also seen that these poles correspond to a 
higher frequency of oscillation than the Langmuir mode. As 
K gets smaller these poles rapidly diminish in width 
compared to the Langmuir pole, move in towards the origin, 
and line up on a line at an angle of 45 degrees. At 
approximately K=0.3 the first of these poles has a frequency 
less than that of the Langmuir mode. AS K grows even 
smaller, more of these poles move in with a lower frequency 


than the Langmuir mode. 


The fact that no resonance exists in E£, for a 
temperature of 100 eV and a density of 4x10!7 cm—3, as shown 
in Fig. 2, even though Fig. 4(c) shows thatw; is close to a 
pole, can now be understood in terms of the form of D(K,w). 
This temperature and density yield, K=0-15. When K is this 
small the higher order poles of the function 1/(K2D) are 
very narrow, and even though they are close to the real 


axis, they do not appreciably alter the shape of the 


function along the real axis. 
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FIG. 5 (a) Contour map of the function 
1K2D (Kw), for K=0.15. Minimum contour 
value is 0.01 and each sucessive contour 
is increased by 0.01. Maximum value is 
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FIG. 5 (b) Contour map of the function 
{K2D(K,w)|], for K=0.50. Minimum contour 
value is 0.04 and each sucessive contour 
is increased by 0.04. Maximum value is 
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FIG. 5(c) Contour map of the function 
{K2D (Kw), for K=2.0. Minimum contour 
value is 1.0 and each sucessive contour 
is increased by 1-0. Maximum value is 
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FIG. 6 (a) Three dimensional view of 
the function 1/({K2D(Kew) 1) as seen from 
30 degrees above the bottom surface. 
Top surface of cube corresponds to a 
value of 100.0. K=0.15 
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the function 1/({K2D(K,w)!) as seen from 
30 degrees above the bottom surface. 
Top surface of cube corresponds to a 


value of 25.0. K=0.50. 


FIG. 6 (c) Three dimensional view of 
the function 1/(|K2D(K,w) 1) as seen fron 
30 degrees above the bottom surface. 
Top surface of cube corresponds to a 
value of 1.0. K=2.0. 
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Figs. 6(a) - 6(c) also show the general dependence of 
E3 ON w3~- From Fig. 6(a) it is seen that a strong resonance 
exists in T/1K2D (Kews) le and thus in Ez, if 3 1s close), to 
the frequency of the Langmuir mode. As K grows larger the 
strength of this resonance becomes smaller. Fig 6(b) shows 
that the resonance is weak when K=0.5. Fig. 6(b) also 
Shows that the value of 1/|K2D(K,w3){ is much larger if w3 
is larger rather than less than the Langmuir frequency. 
Fig- 6(c) shows that almost no resonance exists when K=2.0, 
and that the function, 1/{K2D(K,w3){, is almost independent 
of w3 for a wide range of w3 around the Langmuir frequency. 
Since K is large for high temperatures little importance can 
be placed Bi w3 being close to the frequency of the Langmuir 


mode for high temperatures. 


The general shape of the Langmuir pole can be derived 
analytically for K<<1. The fact that for small kK, the 


Langmuir pole is very close to w3; is used. Now, 


opengl 
Bese w3/k3Vo_ = Te Soe. (2.19) 


Thus for small K, a, _is large. For large % _ 
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Thus, for this case it is found that near the Langmuir pole, 
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KD e K Ue ey fens es 1 

O- O44 (2.20) 
Equation (2.20) has a zero ata 1/ (V2 K) which, by equation 
(2.19), corresponds to yw eo nes Equation (2.20) shows that 
the Langmuir pole Cannes be considered to have a_ simple 
1/(a,_- (1/72K)) dependance when a,. is of the same order of 


Magnitude as K or larger. 


eas COMPARISION WITH INVERSE BREMSSTRAHLUNG 


Figs. 1and 2 show that difference frequency heating 
is relatively insensitive to fluctuations in plasma density 
and temperature. If the powers of the input lasers are 
sufficiently high, a significant enhancement in the rate of 
plasma energy absorption over inverse bremsstrahlung 
absorption can be realized. This enhancement covers a wide 


range of temperature and density. 


The rate of power absorption from the laser through 
inverse bremsstrahlung is given by Johnston and Dawson%, and 


can be written as 


BAyx10ae? wan. wali (Pees 
Fl n n 
2 0 eae (Shee; (2.21) 
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where n is the plasma electron density in cm~3, } is the 
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incident wavelength in cn, Ty is the electron temperature in 
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eV, P; and P2 are incident laser powers in watts/cm2, and 
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The heating rates due to inverse bremsstrahlung can be 


expressed in terms of the dimensionless quantity 


(nee, (BG. te vin ke (2.22) 
0 1/2(P,+P, ) 

where Y2 = 1 sec.watts/(ergs cm2). Table 1 lists values of 
3 Wrpand the ratio W3 /W ip for various values of density and 
temperature. If the power densities of both input lasers 
are equal, (i.e. P|} =P. =P), the ratio in Table 1 multiplied 
by P, yields the ratio of power absorption due to the 
nonlinear process to the power absorption due to inverse 
bremsstrahlung. The inverse of this ratio yields the value 
of P for which W3=Wyp- Table 1 shows that if the power 
densities of the input lasers are both 10!1 watts/cm?, there 
is an enhancement in the plasma power absorption of about 
one order of magnitude for an electron temperature range of 
500 ev to 5000 eV and a density range of 5x10!© cm-3 to 
5x1017 cm. Since the nonlinear power absorption is 
proportional to P?, while inverse bremsstrahlung is 
proportional to P, a plasma power absorption enhancement of 
two orders of magnitude would be experienced if the input 


lasers had power densities of 1012 watts/cm@. 
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The efficiency of the transfer of laser energy to 
plasma energy via this nonlinear method is limited by the 
Manley - Rowe relations. This relation says that if P, is 
the power gained by the wave at frequency w;, and P, is the 
power lost by the wave of the high frequency laser at 
frequency w,, then 

Py ws 

Sor (2.23) 
This relation places an upper bound of #3/w, on the fraction 
of energy of the high frequency laser that can be 
transferred to the plasma through the nonlinear process. It 
ae interesting to note that while the incident beam at 
frequency 0) is heavily damped by this process, the lower 
frequency beam at frequency w, is enhanced. The energy of 
the high frequency beam is transferred to the lower 
frequency beam in such a way that hw, =hw2 thw3. Therefore, 
for maximum efficiency the high frequency laser should be 


more powerful than the low frequency laser. 
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CHAPTER 3 


THE SA. .H5D.) .CODE 


Soil THE HYDRODYNAMIC MODEL 


Fig. 7 shows a schematic diagram for a possible 
experiment on the beat frequency heating of a plasma. 
During such an experiment both the plasma density and plasma 
temperature change at a rate dependent on the heating rate. 
The heating rate in turn is dependent on the plasma density 
and temperature. A prediction of the total energy that can 
be coupled to the plasma through the beat frequency heating 
process can be made with a computer code which simulates the 
MHD behavior of the plasma as it is heated. The MHD code 
used is a modified version of the code described by Burnett 


and Offenberger!°®. 
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FIG. 7 Schematic of a beat frequency heating experiment. 


A configuration of two antiparallel, well focused beams 


mixing in the center of the solenoid is assumed. The laser 
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radiation in the center of the solenoid will cause the gas 
to heat in the center and expand outwards. Since the plasma 
is a good conductor the magnetic field is trapped in the 
gas, and thus expands out with it. The computer routine 
incorporates an exact solution of the appropriate 
magnetohydrodynmic equations to describe the radial dynamics 


of the plasma in the solenoid. 


The following discussion of the routine follows that 
given by Burnett and Offenberger!9°, The only major 
difference between this routine and that used by Burnett and 
Offenberger is that this version separates parallel and 


perpendicular electron temperatures. 


A fully ionized three temperature plasma with a 
parallel electron temperature of 1, a4 perpendicular 
electron temperature of T,, and an isotropic ion temperature 
of T; is assumed. The ions and electrons are assumed to 
have the same local number density, Nos and the same radial 
fluid velocity v. Apart from hydrodynamic terms, energy is 
transported by electron and ion thermal conduction. 
Radiation transport is neglected. The magnetic field is 
assumed to diffuse due to a classical electrical resistivity 
with joule heating being absorbed by the electrons. In 


addition to electron heating due to inverse bremsstrahlung 


absorption, the present version of the routine includes the 


nonlinear heating process discussed in chapter two. The 


power absorbed through the nonlinear process is converted 
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into the component of electron thermal energy parallel to 
the magnetic field. This makes it necessary to calculate 
the perpendicular and parallel electron temperatures 
separately. Parallel and perpendicular thermal energies are 
assumed to be exchanged according to electron-ion and 
electron-electron collision rates. Electrons and ions are 
assumed to exchange energy at the classical equipartition 
rate. An artificial viscosity is introduced to numerically 
handie the Shock wave which is generated by strong 
expansion. The artificial viscosity has the effect of 
spreading the shock front over several mesh points. This 


technique is described by Richmeyer and Morton!!. 


In this code the z-dependence of the variables is 
neglected. The implications of this assumption are 
discussed later. The differential equations that are solved 
numerically to determine the M.H.D. behavior of the plasma 
are given below. 


Continuity or Conservation of Mass 


9 (35-1) 


where r=radius as measured from the center of the solenoid 


and v is the fluid velocity (which is defined to be in the 


radial direction). 
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Conservation of Momentum 


OV OV 9 , 

m.n_(—) + yo - 

Dolan) Jeol 4 Nokg(T, +T,) + P+ Q)= 28 (a) 
o 2 


where Fy = be %% (E2/8T) is the ponderomotive pressure of the 
light and Q is the artificial viscosity term discussed 


above. As given by Richmeyer and Morton1!1 this term is 


Oe=Atmen (as av/ar < 0 


n= 0 av/ar > O 


where A is a constant chosen in such a way as to cause the 
shock front to spread over a distance of about three times 
the width of the spacial net used in the routine. Except at 
a shock front, this term is very close to zero and is 
totally insignificant. For the case of a strong magnetic 
field in the solenoid, shock fronts do not form and this 


term could be dropped from the routine. 
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ej-[e n/(4n) J (By is power from joule heating, 1 is the 
electrical resistivity, Taq is electron-ion thermal 
equalization time, ep is inverse bremsstrahlung power input, 
and Ee is electron collision time. 

Sucessive terms on the right hand side denote work done by 
expansion, thermal conduction losses, joule heating, energy 
transfer to ions, inverse bremsstrahlung energy input, and 
energy transfer from parallel electron energy to the 
perpendicular electron energy. éTy is the difference 
between the electron parallel tenperature and the electron 


perpendicular temperature. (6T, =T,-T, ) It is calculated 


fron, 
2e 
3 vo; ee Nl el ee) 
ep aplcbu) ay Te (3. 3b) 


where ENL is the power input due to the nonlinear heating 


process. 


(3.4) 


The various terms in equation (3-4) denote quantities 
analogous to those defined in equation (3.3a). The first 


term has been nodified by the shock heating. The parameter, 
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Ki is the ion thermal conductivity divided by Boltzmann's 


constant. 


Magnetic Field Transport Anda Diffusion 
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3.2 NUMERICAL INTEGRATION 


The MHD equations are numerically integrated by a 
Lagrangian scheme. The solenoid is divided into 50 to 100 
concentric shells. In a Lagrangian scheme the boundaries of 
the shells move with the fluid velocity and must be adjusted 
with each time increment. The difference equations are 
written in such a way that all quantities except v are taken 
at half-integral Space steps and integral time steps. The 
shell velocity, v, is computed at half-integral time steps 
and integral space steps. The hydrodynamic equations (3.1) 
and (3.2) are treated explicitly and the diffusion equations 
(3-3), (3-4), and (3.5) are treated implicitly. The method 
of solution of the implicit difference equations is 
discussed by Richmeyer and Morton!?. Since the properly 
centered difference equations require knowledge of the 


transport coefficients at half-integral time steps, 


equations (3.3) to (3-5) are solved with an iterative 
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quantities are known, is shown in FiGeeid 


The boundary conditions applied are; 


eas 6 at r=R, 
a9B/ w=0 at c=R, 
v=0 at T=R, (3. 6) 
oT/ or=0 at r=0 
dB/ or=0 at r=0 


These boundary conditions correspond to what would be 
experienced if the solenoid was surrounded by a _ copper 
shell. The plasma temperature at the shell would go to the 
shell temperature. The second boundary condition would 
occur due to the high conductivity of copper and corresponds 
to conservation of Bearer flux. These boundary conditions 
also apply to the case of a larger solenoid, where the laser 


heated region is much smaller than Ros 
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Fige 8 Time and space positions where different 
quantities are known. ‘{' denotes positions and 


times where shell radii are known. ‘'x' denotes 
positions and times where the parameters,T,, Ti ¢ 
T., n., and B are known. ‘v' denotes positions 
ahd times where shell velocities are known. ‘Te 
denotes positions and times where transport 


coeficients are known. 
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3.3 ALTERATIONS TO ROUTINE 


A number of changes have been made to the routine as 
described by Burnett and Offenbergert!®, The program 
efficiency has been increased, and other changes were 
necessary to accommodate the nonlinear heating process. A 
small error was found in the calculation of the joule 
heating term of the original routine, but its correction 
made an insignificant difference to the results. Also 
changes made in the calculation of the transport 
coefficients made little difference to the final results. 
The contribution to heating from inverse bremsstrahlung was 


changed to that given by Johnston and Dawson?. 


One major change already discussed was the inclusion of 
(6T,) Seperate from the electron temperature, which was made 
necessary by the nonlinear heating process. The factor, Tor 


which is the electron collision time is calculated fron 


Nags Ne piel (3.7) 
Te Tej Tee 
| where Da is the electron-ion collision time and eek the 


electron-electron collision time. Tog is given in (5.26) of 


Spitzer!2, and Ta] Can be calculated from (5.22) of Spitzer. 


Another major change, which greatly increased the 


efficiency of the routine, was the inclusion of a subroutine 


which made the program adjust its time step size to the 


maximum size allowed by the plasma conditions at that time 
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step. 


This was accomplished by choosing the time step, At, 
such that the following conditions are satisfied: 

1- At must be less than a specified fraction i debe 
the maximum time step allowed by the stability 
condition. (The stability condition is discussed 
below.) 

2. The radius of any shell boundary cannot change by 
more than ten percent in a single time step. 

Se eeuhe shell thickness of any shell cannot change by 
more than ten percent in a single time step. 

4. The electron or ion temperature cannot change by 
more than ten percent in a single time step. 

5. The size of a time step cannot increase by more 


than fifty percent in a single time step. . 


The size of the time step used is usually limited by 
the first condition above. The stability condition used was 
one derived by Richmeyer and Morton!! and can be written as: 

2 oe ale 

At < A, Mintar,(C A (3-8) 
where C is the adiabatic sound speed and Cp is the Alfven 
velocity. Essentially all this condition says is that the 
time step must be smaller than the time necessary for a 


pressure wave to travel a distance equal to the smallest 


Space step. 
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The stability condition was derived for a case in which 
no input power was given up to the plasma. Equation (3.8) 
is found to be quite adequate when the laser power input is 
low if A, is set to 0.8. When the laser power input is 
large (i.e. an input laser power density of 1012 watts/cm2) 
or when there is good nonlinear energy coupling between the 
laser and the plasma, it is found that AY must be reduced to 
0.5. For a larger AY the routine will generate oscillations 
near the center of the solenoid, with a time period equal to 
the time step size and a space period equal to the space 
step size. For economy it is desirable to make A, as large 
as possible. It is found experimentally that the results 
are unaffected by the choice of A, as long as it is small 


enough to prevent the oscillations previously mentioned. 


Finally the boundary condition for B at r=R was 
changed from B=0 at r=R) as in the original version of the 
routine to 9B/ar =0 at r=R)- This is easily derived as 
being the boundary condition for a cylinder surrounded by a 


copper sheet with the use of Ampere's Law. 


This last change implies conservation of magnetic flux 
in the routine. It allows the use of conservation of energy 
and conservation of magnetic flux as two checks on the 
numerical procedure. It is found that in a 40 ns run both 
of these quantities are conserved to 5 Significant digits. 
For longer runs the temperature starts to rise at the 


boundaries and energy can be lost from the plasma. 
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CHAPTER 4 


RESTRICTIONS AND ASSUMPTIONS 


4.1 SOLENOID LENGTH 


For the one-dimensional M.H.D. routine to give a 
reasonable approximation to reality, it must be assumed that 
the solenoid's length is much greater than its radius, and 
that the heating is reasonably constant over its length. In 
the M.H.D. simulations made for this thesis the radius of 
the solenoid was assumed to be 0.5 cm, so if its length is 
about 5 cm or greater the first of the above conditions will 
be satisfied. The model may be valid for a shorter solenoid 
if calculations are made for a short enough time to prevent 
Significant disturbances to reach the outer radius of 0.5 


CMe 


The heating in the solenoid is greatest close to where 
the high frequency laser enters the plasma. Here both 
lasers would have maximum intensity; the first one because 
it hasn't had time to damp, and the second because its power 
grows as it propagates through the plasma. The power lost 
by the first beam is approximately ten times as great as the 
power absorbed by the plasma, and the power gained by the 
second beam is approximately nine times aS great as the 
power absorbed by the plasma. Using this information the 


length over which the product P) P) is reasonably constant 


can be calculated with the use of Table 1 of chapter two, 
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for specific values of Pie Poe Ty, «60cand 6 le The power 
O 
intensity of the first laser P, must satisfy the equation 


oP, 
ax, P(x) PDAx) (4.1) 


where a=10—6 Wen. The constant a is equal to ten times the 
power absorbed by the plasma in watts/cm3 divided by the 
product of the laser intensities in watts/cm@. Assuming the 
laser beams are insignificantly damped by inverse 
bremsstrahlung absorption, P(x) can be written in terms of 


P,(x) as 
Paix) "Pa (x) + Po(=).- P, (a) (4. 2) 


where P.() is the power of the second laser beam before it 
enters the plasma and P, (~) is the power of the first laser 
beam. after it has passed through the plasma. Using 


equation (4.2) the solution to equation (4.1) is 


| b P,(o) exp[-bx] 


4.3 
Paix) a P;(o) (1 - expl-bx])+b art 
where =a(P>(~)-P,(~)) and P, (0) is the power of the first 


laser before it enters the plasma- With this solution for 
P, (x), Py (x) can be calculated from equation (4.2). 


Approximate values of q for different temperatures and 


densities can be found with Table 1 or with Figs. 1 and 2. 


Equations (4.2) and (4.3) are used in chapter 5 to find the 
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distance over which P,P, is reasonably constant for the 


conditions experienced in the MHD simulation. 


4.2 STIMULATED SCATTERING 


When an intense electromagnetic wave is present in a 
plasma, two new waves can be parametrically excited by the 
incident or pump wave. If one of the excited modes is 
electromagnetic it can escape from the plasma and appear as 
stimulated scattering of the incident wave. The process of 
stimulated scattering is very closely related to the 
nonlinear heating process studied in chapter two. In the 
process of stimulated scattering, the incident wave of 
' frequency Wo and wavenumber K | interacts with a propagating 
density perturbation (w,k) associated with an electrostatic 
wave. This results in an electromagnetic wave of frequency 
ow and wave vector k -k being scattered. This scattered 
wave can then interact with the incident wave in much the 
same way as the waves from the two lasers described in 
chapter two. This interaction then amplifies the original 
density perturbation. An unstable situation can result if 
‘the energy transfer to the scattered electromagnetic wave 
and the electrostatic wave exceeds their natural damping 
rate. Some energy is deposited in the plasma via the 


electrostatic wave, but the amount is limited by the 


Manley - Rowe relations. 


Theoretical studies of this process, assuming all waves 
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are coherent have been made by Drake et al.13, A similar 


study, but assuming incoherent waves has been made by 


Tsytovich!4, 


Theory shows that the incident wave intensity has to 
exceed a threshold value before significant scattering will 
occur. After the threshold intensity has been exceeded, the 
scattered wave amplitude grows exponentially according to 
e’t, where y is known as the growth rate and t is the time 
over which the wave can grow. The time t is limited to 


t=2/c, where 2 is the interaction length and c is the speed 


ofelight. 


Stimulated Brillouin scattering (the scattering of 
light from ion-acoustic waves) could possibly present a 
problem in a solenoid such as the one described in this 
thesis. Experiments have been made on stimulated Brillouin 
scattering by Offenberger, Cervenan, Yam, and Pasternak!®%. 
The Manley -— Rowe relations tell us that very little energy 
can get into the plasma through stimulated Brillouin 
scattering since the ion-acoustic wave frequency is small 
for the values of wave vector, k, which occur in such a 
process. For example, if the electron temperature is 100 eV 
and the plasma density is 4x10!7 cm~* a maximum of the 
fraction 6-4x10-% of the incident beam's power can be 


transferred to the plasma through this process. 


The theory from Drake et al.13, and the experiment by 


Offenberger et 1.15 both indicate that for power 
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intensities. of 1014! to 1012 watts/cm2, and for temperatures 
and densities assumed in this thesis, growth rates for 
stimulated Brillouin scattering of about 1011 to 1012 sec-1 
can be expected. Such growth rates over several centimeters 
would ensure that almost all the light would be scattered. 
However the interaction length, 2, over which the scattered 
beam can grow will not be larger than the coherence length 
of the ion-acoustic wave from which it scatters. Ina 
rapidly heated plasma a certain amount of turbulence would 
be expected. This would limit the coherence length of the 
ion-acoustic wave. Experimental evidence given by 
Offenberger et al.15 indicates that the coherence length of 
the ion-acoustic wave may be sufficiently small, so that 


only a small amount of light will be scattered. 


The process of scattering from electron plasma waves, 
or stimulated Raman scattering, is also very closely related 
to the heating process of chapter two. In this process the 
scattered wave would propagate antiparallel to the incident 
wave, and its frequency would match the electron plasma 
frequency. The Manley - Rowe relations would impose 
approximately the same maximum energy transfer condition 
that was imposed in chapter two. This process would appear 
also as a method of enhancing heating. The usefulness of 
this process to enhance heating is limited by the fact that 
the growth rate of stimulated Raman scattering is smaller 
than the growth rate of stimulated Brillouin scattering for 


the plasma conditions assumed in this thesis, and by short 
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natural coherence lengths of plasma waves. 


4.3 DISTORTION OF THE VELOCITY DISTRIBUTION FUNCTION 


The equations which describe the nonlinear heating were 
derived assuming a Maxwellian velocity distribution 
function. It might be expected, however, that the continual 
Landau damping of the driven wave would cause the 
distribution function to flatten near the phase velocity of 
the difference frequency wave, w,/k;. This flattening of 
the distribution function would cause a decrease in the 


Landau damping. 


It turns out that the distribution function remains 
much closer to Maxwellian than would be predicted as a 
result of collisions. Theory shows that collective 
electrostatic interactions have an effect similar to 
collisions on the distribution function, and that a thermal 
velocity distribution can be maintained, even in the absence 
of collisions. Bernstein and Trehant® show how a limit is 
placed on the magnitude of plasma oscillations in a cold 
plasma. At this limit a process Similar to the breaking of 
water waves takes place. This process causes ordered motion 
to go over into unordered thermal motion. One dimensional 
many body computer models (see Bernstein and Trehan!®&) allow 
this process to be investigated in warm plasmas. It is Seen 
that ain this case, energy is fed even faster into thermal 


motion from the organized wave energy. A discussion of this 
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process can also be found in chapter 7 of Schmidt17 


Computer models by Dawson!8, and Dawson and Shanny?9, 
show that the distribution function will maintain a 
Maxwellian form in the presence of Landau damping much 
better than would be predicted on the basis of collisions 
alone, but af the amplitude of the plasma wave is 
sufficiently high, the damping rate will deviate from that 


predicted by linear Landau theory. 


Most of the theoretical work done in this area has been 
done as an investigation of the time development of finite 
initial value amplitude electron waves, rather than the 
steady state solution of a driven wave as is the case at 
hand. Three approaches have been used to investigate the 
time development of these electron waves. The results of 
these have been summarized in an experimental paper by 
Franklin, Hamberger, and Smith2°, The first approach is 
that of the one dimensional computer model and has already 


been discussed. 


The second method, used by Armstrong?!, is to solve 
Simultaneously Vlasov's and Poissons's equations using 
numerical techniques. The results, summarized by Franklin 
etemeat.-9, Sacre that? the damping rate, y, , is correctly 
predicted by linear Landau theory if the amplitude of the 
condition can be expressed 


plasma wave is. small. The 
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: << el 
cel (4.4) 


¢ is the potential of the electrostatic wave. When () 
increases, violating equation (4.4), the damping rate 
decreases. When eo/Kk ,1™1 the damping rate is no longer 
exponential, and the instantaneous rate exceeds the value 
predicted by linear Landau theory. This behavior can be 
associated with the "wave breaking" process in a cold 


plasma. 


The third method separates the electron distribution 
function into a resonant part and a non-resonant part, and 
solves exactly the Vlasov equation in the resonant region. 


This technique was used by O'Neil?2, and later by Bailey and 


Denavit23. Their results showed a decrease in Landau 
damping if 
Ke 

SE erage ee (4.5) 
B 

where uw, =k; (e9/m, ) is the frequency of oscillation of the 


trapped electrons. kK3 and » are the wave vector and 
potential respectively of the electrostatic wave. Sugihara 
and Kamimura?2* computed self-consistent equilibrium 
solutions to the initial value problem for a range of 


O<q Kio They found the same result as O'Neil for q~ 0 and 


linear Landau damping for q * ®- 
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An experiment to test the damping behavior of plasma 
waves for various values of eo/k, Ty and q was done by 
Franklin, Hamberger, and Smith2°,. They did the experiment 
on a sodium column, with an electron temperature of 2500 °K 
and a density of about 2x107 cm-3. The experimental results 


indicate that linear Landau damping could be expected for 


poe < exio? (4.6) 


It should be noted here that this condition was derived and 


tested only for the case where K is of the order of unity. 


When K is of order unity it is assumed that relation 
(4.6) is an approximate condition which must be imposed on 
the difference frequency wave before the equations for the 
nonlinear heating are to be mecurate: For the specific case 
of difference frequency heating with lasers of wavelengths 


of 9.6 and 10.6 microns condition (4.6) becomes 


ete SOA Gey IO (4.7) 
Values of E3; for different values of nos T, and laser 
powers can be found from Figs. 1 and 2 of chapter two. When 
using this condition to find restrictions that must be 
placed on the routine, it must be remembered that by the 
time the laser power, and thus £3, reach a maximum, the 


temperature of the electrons is high. It should also be 
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remembered that this condition was derived by neglecting the 
role of collisions in thermalization. at low temperatures 
collisions can play an important role in maintaining the 
Maxwellian distribution. If the electron temperature is 
sufficiently high to make K of the order unity, relations 
(4.6) and (4.7) can be used to check the validity of the 
results. For example, at a density of ee SLES an 
electron temperature of TT, =1000 eV, and input laser 
intensities of 2.5x101!1 watts/cm2 it is found that E, =30 
statvolts/cn. For hic) temperature and density, 
K=kAp =0.46. Relation (4.7) then yeilds eo/k, Ty =7- 2x10-%. 
This implies that the driven wave should be experiencing 
damping as predicted by linear Landau theory. If the 
electron temperature is lowered significantly from this 
value, K will decrease, and relation (4.7) cannot be used to 


test the validity of the theory. 


The effect of E3 exceeding the restrictions placed on 
it by relation (4.6) can be predicted qualitatively, with 
the aid of Figures 3 and 4 of chapter two. When E3 gets 
large the Landau damping will decrease and the poles 
associated with D(K,w) will move towards the real axis. If 
w3 is very close to the frequency of the Langmuir mode, a 
stronger Saeley resonance could be expected leading to a 
greater heating rate than that predicted in chapter two. 
If, on the other hand, w3 is not close to the frequency of 
the Langmuir mode, the motion of the poles towards the real 


axis would have little effect on the value of F;, but the 
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decreased damping would lead to a reduction in power 


absorption. 
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CHAPTER 5 


RESULTS 


The M.H.D. code discussed in chapter 3 has been used to 
compare the heating of a plasma when two lasers are mixed 
with a difference frequency near that of the Langmuir mode, 
with the heating that is obtained when inverse 


bremsstrahlung is the only energy absorption mechanisn. 


The laser energy pulse assumed in the routine is 40 
nsec long and of triangular shape. The pulse starts with 
zero power, reaches peak power at 20 nsec, and falls off to 
zero again at 40 nsec. The radial intensity distribution of. 
the focused laser pulse is assumed to be GausSian. The 


intensity as a function of radius, rt, can be expressed as, 
veil ji 
I(r) =P /(2no )exp[-r*/20°] (5.1) 


where P, is the laser power, and o is the radius of the 
beam. o has been set of 0.045 cm. For all the runs 


reported in this thesis both lasers were assumed to be of 


equal power. 


The spatial mesh in each case consisted of 50 shells. 

Initially the inner shell had a radius of 3x10-3 cm. The 
: ‘ = n 5 

thickness of the n'th shell is given by 3x10-3xA where A is 


chosen in such a way as to make the radius of the outer 


shell 0.5 cn. 
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Initially, a fully ionized hydrogen plasma with 
electron and ion temperatures of 1 eV is assumed. The 
plasma has a uniform electron and proton density of 7x10!7 
cm—3. The magnetic field is longitudinal and uniformly 
distributed throughout the solenoid. Several runs were made 


With different initial magnetic fields and different laser 


powers. 


Pigurecs 9 towlzearesplotssof Ce, Ue, mele 


oe B, and n 


9) 
versus radius at a time of 40 nsec. In Figs. 9(a), 10(a), 
11(a), and 12(a) energy is assumed to be absorbed by inverse 
bremsstrahlung only, while in Figs. 9(b), 10(b), 11(b), and 
12(b) it is assumed that energy is absorbed by both inverse 
bremsstrahlung and the nonlinear process. In all figures 
the plotted value of magnetic field is normalized to the 


initial plasma density. The plotted value is given by 


B/B, x initial plasma density. 
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(b) Nonlinear + inverse bremsstrahlung. 


Fig. 9 Plasma parameters Ty and T eases 
curve), n, (solid curve), and B normalized : 
initial density (dashed curve) as a function 0 
radius at 40 ns. Initial Rech 4 are 

0 = ea aser 
n =7x1017 cm 3, Beeo0 kG, an p 


pdwer=1.25 GW. 
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Fig. 10 Plasma parameters T, and T _ (dotted 
curve), Mm (solid curve), and B normalized to 
initial density (dashed curve) as a function of 
radius, at 40 NSe Initial conditions are 
n.=7x1017 cm-3, Boel e0 kG, and peak laser 


power=1.25 GW. 
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Plasma parameters T, andi (dotted 
n (solid curve), and 8B normalized to 


density (dashed curve) as a function of 


at 40 NS. Initial conditions are 


Dope 1007 cm-3, 3B. =100_ kG, and peak laser 
power=12.5 GW. 
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Ins Figs. 9 (a). and 9(b) the following parameters were 
used; initial magnetic field, Bi =60 kG.; initial plasma 
density, Det XAOte cm—3; peak laser power, P)=1.25 GW. 
Total laser energy is 25 joules. Comparing these two 
figures, it is seen that the temperature of the electrons in 
the center of the solenoid is about 330 eV after the 
nonlinear heating mechanism is used, and about 200 eV after 
only inverse bremsstrahlung absorption is used. For this 
run the parallel and perpendicular electron temperatures 
remained within ten eV of each other throughout the 
Simulation, and were essentially equal at 40 ons. It was 
found that in the center shell the ratio, W3/Wipe varied 
from 5 to 8 during the time period of 18 ns to 35 ns, where 
W3 is the heating rate due to the nonlinear process and Wp 


is the heating rate due to inverse bremsstrahlung 


absorption. 


Figures 10(a) and 10(b) come from a simulation where 
the parameters are all the same as those in Fig. 9 except 
for the magnetic field which has been increased to 100 kG 
for this run. Increasing the magnetic field makes it more 
difficult for the plasma to expand outwards during the 
simulation. This results in the density staying higher in 
the center of the solenoid, which increases the coupling of 
laser energy to the plasma through both inverse 
bremsstrahlung and the nonlinear processes. This is shown 


in Fig. 10 where the electron temperatures and plasma 


densities are higher than those shown in Fig. 9. Again, the 
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electron temperature in the center of the solenoid is 


Significantly higher for the case where nonlinear heating is 


assumed. 


Bethe and Vlases25 have suggested that one of the 
requirements to heat a plasma to thermonuclear temperatures 
in a long solenoid is a magnetic field of 300 to 500 kG. 
Thus, it is of interest to make a simulation of a laser 
heated solenoid with fields of this strength. Figs. 11(a) 
and 11(b) show the results of such a simulation, at a time 
of 40 ns. Again, the laser power is the same as in previous 
runs; total laser energy is 25 joules, peak laser power is 
1.25 GW. These plots show that the density stays almost 
constant throughout the simulation. The nonlinear process 
is able to heat the electrons in the center of the column to 
about 1502.5 times as high a temperature as inverse 


bremsstrahlung absorption alone. 


As explained in chapter two, the power absorption due 
to the nonlinear process is proportional to the product of 
the powers of the two lasers, while the power absorption due 
to inverse bremsstrahlung absorption is proportional to the 
sum of the powers. Thus W3/Wrp increases as laser powers 
increase. The effect of increasing the “iweer, power was 
investigated by making a simulation where the lasers had a 
peak power of 12.5 GW and a total energy of 250 joules. The 
initial magnetic field was set to 100 kG. The results of 


this simulation are shown in Figs. 12(a) and 12(b). The 
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final electron temperature for the case of nonlinear heating 
is more than three times the electron temperature with only 
inverse bremsstrahlung absorption. Fige 12(b) shows that 
the parallel and perpendicular electron temperatures have 
not yet equalized at 40 ns. It is interesting to note the 
sharp density minimum in the center of the solenoid 
(Fig. 12(b)) caused by the rapid nonlinear heating. In the 
center shell the ratio, W3/W pe was between 8 and 40 for the 


period of 5ns to 35 ns. 


So far no mention has been made of the resulting ion 
temperatures in these simulations. Fig. 13 shows a plot of 
the ion temperatures for all the simulations made. These 
plots show that in the more rapidly heated plasmas, 
decoupling between electrons and ions is greater than in the 
plasmas that are heated more slowly. These plots also show 
an advantage in having a large magnetic field. The large 
magnetic fields maintain a higher plasma density in the 
center of the solenoid, and this makes the electron-ion 
equipartition time smaller. Spitzert?, shows that the 
equipartition time is approximately proportional to ei diee 
The very rapid heating experienced in the simulation with 
the 250 joule laser and the nonlinear heating resulted in 
almost complete decoupling between electrons and ions in the 
center of the solenoid. This is seen as a steep minimun in 


ion temperature in the center of the solenoid. 
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The value of 2 9/k, Ty has been calculated at several 
times in the simulations. In general, the value of e¢/k, T,, 
is so close to the critical value defined in relation (4.7) 
that it is impossible to predict the accuracy of the damping 
rate of the beat frequency wave, as given by linear Landau 


theory, under these conditions. 


For example, in the simulation where the laser energy 
is 25 joules and the initial magnetic field is 100 kG, 
K=0.17 and e6/kpT, =4.3x10-3 at 15 ns; at 35 ns, K=0.3 and 
ea =5.8x10-*. Similarly for the simulation with the 
higher powered lasers (laser energy of 250 joules, and 
initial magnetic field of 100 kG) K=0.68 and 
eo/kpTy =3.4x10-3 at 15ns; and K=1.3, and ed /Kkp Ty =5 25x 10=4 
at 30 ns. For both of these simulations, the values of 
eo/kpT, found at 15 ns is larger than the critical value 
defined in relation (4.7). However, Since relation (4.7) is 
only an approximate condition it is difficult to estimate 
the accuracy of the nonlinear heating equations given in 
chapter two under these conditions. Relation (4.7), 
however, does predict that the calculated nonlinear heating 


rate should be accurate under the conditions experienced at 


3oense 


The total energy absorbed by the plasma can be 
calculated by adding up the total thermal, kinetic, and 


magnetic energy in the plasma initially, and after 40 ns. 


The absorbed energy is expressed in ergs/cm length of 
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solenoid. The increase in total absorbed energy due to the 
nonlinear process is Surprisingly small. For the simulation 
with 25 joule lasers and 60 kG Magnetic fields the total 
power absorbed was 1229x107 ergs/cm for inverse 
bremsstrahlung absorption alone, and 1.34x107 ergs/fcm when 
the nonlinear process was used. This is an increase in 
total energy by a factor of only 1.04. This is partially 
caused by the fact that power absorption from both 
mechanisms decreases as the temperature rises and density 
falls. Another factor to be considered is that, because the 
nonlinear power absorption is proportional to the product of 
the laser powers, the heating takes place over a smaller 
cross sectional area of the solenoid than does inverse 
bremsstrahlung heating. It must also be remembered that the 
thermal energy is proportional to the product of plasma 
density and temperature. Thus, even though the temperature 
is much greater in the center of the solenoid for the case 
in which nonlinear heating is used, the total thermal energy 
in the center of the solenoid is not that much greater, 


since the plasma density is lower. 


In the simulations where the initial nagnetic field is 
100 kG and laser energies are 25 joules, 1.59x107 ergs/cm 
a absorbed when only inverse bremsstrahlung absorption is 
used, and 1.69x107 ergs/fcm are absorbed when nonlinear 
heating is used. For this case the nonlinear heating 
results in a total energy absorption enhancement of a factor 


of only 1.06. For the Simulations with a 300 kG field, 
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1.89x107 ergs/cn are absorbed for the case of only 
accounting for inverse bremsstrahlung absorption and 
2.-06x107 ergs /cm for the case including nonlinear heating. 
This represents a total energy absorption enhancement of a 


factor of 1.09. 


The simulation with 250 joule lasers and an initial 
field of 100 kG shows that nonlinear heating results in 1.5 
times aS much energy being absorbed by the plasma. With 
only inverse bremsstrahlung absorption, 3.2 8x107 ergs/cm are 
absorbed, while 4.89x107 ergs/cm are absorbed when the 
nonlinear process is used. Thus, if the laser powers are 
sufficiently high the nonlinear process can produce a 


significant enhancement in total energy absorption. 


These simulations are most accurate for a solenoid of 
about 3 cm in length. If it is much shorter than this, the 
one-dimensional approximation is no longer valid. If it is 
much longer than this, there would be times in the 
Simulations when the high frequency laser beam would be 
depleted at one end of the solenoid. Even for the case of a 
3 cm solenoid there are times that the nonlinear process 
deposits much more energy at one end of the solenoid than 
the other. For example, in the run with the 25 joule lasers 
and a magnetic field of 100 kG, the high frequency laser is 
depleted to about u0% of its initial value in 3 cm, at a 
time of 16 nse The low frequency beam grows to about 150% 


of its initial value. Consequently, at one end of the 


nr 
yw 4 | 20 sitet tena sna ‘eexeen sae 


f a » 


ws! - ; 
rs ZA T 14aek sive Ors iv g 2 
7 7 : t 2 vaks aan a 145 70 


ok «ah 9th val dind tetas Pace inf9 Ss 9 ‘oor te epi 


a 4 


9 Ria ed? Fil paar emis Sei Apa 2s 
4 peifasocto vo io caas ree barr 


: we 14 Jw. J ' anaes sanurt P| white 


Meni | Hadad wh? ¥) a uw? «year A senor — sa 
—— runes sen: } ae | ‘ae efsastels 


Yoseds Lejos a tatiana teok 
a 
- 
jo DiGasivs r S3650C0s -Fawd oF RRR vende 


= an ) Heald weltpge, ene i 22 FF daguin coe oe >c 


sav 9AeWSh Se os. nel adbscuapa angie: the 


;. 2@ef" (of “Digi ‘seada! -caddd ae age SS: 
7 


3 5 j 7. oe ¥O0 rs e349 doid ay smile «PRG. He sivel 
: — ae 
» toners) ind perry niGhe [aFe8d? Sok ne Seatalh 
) a a - * 
‘qheqad)- teppliach aif? yess eseis ote Seed & nolo 5 
ep aire canis os Some Stella 


« 


sof? Btereicnn 64° 20 bhaleno 2s cart . ae coal ba 1) 
‘ ' 


. ; 
epeud) @f§our $9. etr atiy oes Ail mt .% en at _ 
Gas aca 
eer ete 


hb ss mak Anmainaiy 3) ee i 
4 A? + : 
© 28 4? hal ‘Geter ime Int ad a rs 


7 r 
: _ 
‘ , i Di@ } > 
0 Fs #4 os 


— 


—— 
~— ws a byt i 


oman oniey 
mares Sem! 
_ 


67 


solenoid the actual heating rate is 150% of the rate used in 
the simulation, while at the other end it is only 40% of the 
rate used in the simulation. At this time in the 
Simulation, the electron temperature in the central shell of 
the solenoid is about 200 ev. The parallel thermal 
conductivity (Calculated from equations (5.47) and (5.48) of 
Spitzer!2) is 4.2x107 watts/cm eV. The calculated value of 
power being deposited in the plasma is about 2.3x109 
watts/cm3. Thus, it may be expected that one end of the 
solenoid will start to heat faster than the rate calculated 
in the simulation. But, if the temperature increases 
Significantly, the nonlinear coupling between the laser and 
plasma decreases, allowing the laser to penetrate further 
into the plasma. The thermal conductivity is then 
sufficiently large to keep the temperature of the solenoid 
reasonably uniform along its length. At later times in the 
Simulation there is no problem in maintaining a constant 
temperature along the length of the solenoid since the 
thermal conductivity rises according to T 0/2 and the high 


frequency laser penetrates the plasma more easily as the 


temperature rises. 


For the simulation with higher powered lasers, the 
parallel thermal conductivity is much higher due to the 
higher electron temperature. AS a result a fairly uniform 
temperature distribution is maintained along the length of 
the solenoid, even when the first laser beam is depleted 


significantly as it propagates through the plasma. 
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CHAPTER 6 


CONCLUSION 


The analysis in this thesis indicates that the mixing 
of two laser beams to generate a beat frequency harmonic 
near the electron plasma frequency can provide a useful 
technique of heating the electrons in a plasma. It has been 
Shown that this process is insensitive to changes in density 


and temperature. 


At low temperatures a serious limitation in the theory 
arises due to the difficulty in estimating the accuracy of 
the damping rate of the beat frequency wave as predicted by 
linear Landau theory. However, at electron temperatures of 
300 eV and higher, theory predicts that the damping rate 
will be very close to that predicted by linear Landau 


theory. 


The results of the M.H.D. simulations show that this 
nonlinear process can greatly enhance the heating of 
electrons in a short experimental solenoid. In particular, 
sip the laser beam intensity reaches values of 101? 
watts/cm2, a Significant fraction of the laser energy can be 


deposited in a short solenoid. This is demonstrated in the 


Simulation where two 250 joule laser pulses raise the 


electron temperature in the center of the solenoid to over 
2000 eV. When one considers that the absorption length of 


; 3/2 
inverse bremsstrahlung absorption has a T{° /n, dependence, 


it is seen that it would take a much higher powered beam to 
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achieve this temperature if inverse bremsstrahlung were the 


only energy absorption mechanisn. 


The results of this thesis suggest further work should 
be done in the investigation of the use of beat frequency 
heating of a plasma in a long solenoid. An example of how 
beat frequency heating may lead to enhanced absorption in a 
leng solenoid is given below. A laser fires a high powered 
beam of radiation with a wavelength of 9.6 microns down the 
center of a plasma filled solenoid. Just before this beam 
reaches the far end of the solenoid, another laser fires a 
short, but high intensity pulse of 10.6 micron radiation 
down the solenoid antiparallel to the first beam. When this 

pulse meets the beam of 9.6 micron radiation it will 
interact with it. The 9.6 micron beam will be depleted, and 
the 10.6 micron pulse will grow in amplitude. Thus, about 
90% of the energy of the 9.6 micron beam is converted into a 
pulse of 10.6 micron radiation travelling in the opposite 
direction to the original beam and 10% of the energy is 
naneaieed in the thermal energy of the plasma electrons. 
Turning the beam around like this makes the laser radiation 
travel through a greater length of plasma so that more 
energy will be absorbed by inverse bremsstrahlung 
absorption. However, if all the energy of the 9.6 micron 
beam becomes compressed in the first few cm of the 10.6 
microm pulse as this theory would suggest, the intensity of 


this pulse would grow to such a size that other nonlinear 


processes may become dominant. 
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APPENDIX 1 


COMPUTER ROUTINES 


This routine and associated Subroutines constitute 
the main M.H.D. code, 


aAgaANANAaAAaAA 


INTEGER P,PRTPR 
COMMON P0,BZ,R8,2,T,TF 

FOR ITAPE=1 OUTPUT IS ON TAPE. SPECIFY FILE NAME IN CONTROL CARDS. 
NPRT IS THE RATIO OF (TAPE OUTPUT) /(PRINT OUTPUT). 

NTEST (=0 FOR ORIGINAL START AT TIME C) (=1 FOR RESTART WITHOUT 
CHANGING SHELL SPACING) (=2 FOR RESTART REORGANIZING DATA) 
NREGeIS THE RECORD NUMBER TO BE READ IN FOR RESTART OFF UNIT 4, 
NCYCH IS ZERO FOR INVERSE BRENSTRAULUNG ONLY, IS 1 FOR 
NON-LINEAR HEATING. 

W=3.0E-03 

IT=3 

P=50 

VI=0.9 : E Be 
BZ=3.0E+05 

P0=10.0 

DT=0.25E-10 

TF=40,0E-09 

PRTFR=1 

ITAPE=1 

NPRT=20 

NTEST=1 

NREC=1 

NCYCH=1 

CALL ALAS2 (VI,W,DT,IT,PRTFR,ITAPE, NPRT,NTEST, NREC, NCYCH) 

STOP 

END 
SUBROUTINE ALAS2(V1I,W,DT,IT, PRTFR,ITAPE, NPRT, NTEST ,NREC,NCYCH) 

i EGER P,O PRTFR 

Sites,” £0 (101) , U0 (101) , PEO (100) -D70 (100) ,BC (100) ,¥O (100) 
1,81 (101) ,01 (101) , TEV (109) »B1 (100) ,DT1 (100) ,N1 (100), DN (190), 
2 Ak (109) BB (100) ,CC (100) ,DD(100) »E(100) ,F (100), 

PA(100),DR(100) ,EE,FF, ; ) 
celacet. ani (100) , bb (400) ,cB 190) ,005 100) « 
+) , BBI (100) ,CCI (190), DDI (100) ,AAE (100) , BBE (100) , , 
DDE (160) sD’ 1 (100) ,DR2{100) , DR3 (100) ,R4 (OC) ,D41 (100) -b42 (100), 
7043 (100) , DX4 (100) ,DX5(100) , DX6 (100) DK 7 (100) » PE (100) » ARC 0), 
L1,T10 (100) ,HJUL (100) 

OREAL  W2 (100) sTE2 (100) , P12 (100) -DT2(100) ,B2 (100) «LIN (100) ,WEAT (1 
4 00) , HRATIO (100) /10C¥C .0/yRC (100) ,KE2 (100) ,KI2 (100) ,TEQ2 (109) 
2 KR2 (100) gETHO (100) , EKEO (100) , BBO (100) ,PB (100) , PKE (100) , PABS (100) 
3° PTH (100) ,MI,KB ,TPO(100) , TPF (100) ,TP2 (100) , £3224 (100) , POWER (100) 


REA AOR oA) © MD r 
ake BF 
vi wr vipe 
ae back 
i gs pd ; 


a 


cat 2 nape 

rot pied 3 7 ef, y b¥ oo 
banter ae 7 te ey 
apna SOREN. bie 


east cleus ta 


C10) OGG A) Quasi Qa fay) Cara OOK a i an a 


4,W3N1(100) , TPARD (100) , W3M (100) 

COMMON PO,5Z,R,P,T,TE/COMTHO/TEO, DTO, NO, BO/COMTH1/TE1,DT1,N1,B1 
1/COMHYC/U0,RO /COMHY1/U1,R1/COMTH2/TE2 ,N2, B2 /COMLAS/LIN, HEAT 
2/CONEN/DN,ETHO,EKEO, EKO, PABS, PB, PTH, PKE, PABST, PTHT, PKET, PBT 
3/PARAL/TPO,TPF,TP2,£322, W341, TPARD,W3M . 

SUBROUTINE ALAS WILL SOLVE THE TWO FLUID MAGNETOHY DRODYNAMIC 

EQUATIONS IN CYLINDRICAL SYMMETRY ASSUNING QUASINEUTRALITY BY 

MEANS OF A LAGRANGIAN DIFFERENCE SCHEME 

THE EQUATIONS FOR THE MAGNETIC FIELD AND ELECTRON AND ION 

TEMPERATURES ARE A NONLINEAR COUPLED SET AND ARE SOLVED 

BY AN ITERATIVE PROCEDURE WITH THE TRANSPORT COEFFICIENTS 

EVALUATED AT TIME T+,.5DT COMPUTED ANEW AT THE START OF EACH 

ITERATION 

bZ DENOTES AN EXTERNAL MAGNETIC FIELD 

PO DENOTES THE FILLING PRESSURE(HYDROGEN MOLECULES) 

VI DENOTES A VON-NEUMAN VISCOSITY PARAMETER (SHOCK THICKNESS) 

P (<100) DENOTES THE NUMBER OF CYLINDRICAL SHELLS USED iN THE INT. 

R REPRESENTS THE OUTER LIMIT OF THE SPATIAL MESH 

W REPRESENTS THE THICKNESS OF THE CENTRAL CORE 

DT REPRESENTS THE TIME INCREMENT (HELD CONSTANT IN ALAS2 ) 

TF REPRESENTS THE TOTAL DURATION OF TNE INTEGRATION 

IT DENOTES THE NUMBER OF ITERATIONS IN THE LOOP FOR TE1,D0T1,B1. 

PRTFR = PRINTING FREQUENCY, 

THE FOLLOWING PRECAUTIONS MUST BE OBSERVED WITH ALAS2: 

(1) IF DT IS CHOSEN LARGER THAN THE STABILITY LIMIT AN INSTABILITY 

WILL RESULT IN NEGATIVE DENSITIES AND AN ARITHMETIC ERROR 

WILL TERMINATE EXECUTION. IN GENERAL IT iS SAFE TO WORK WITH 

DIR ERLOVONE-HALE OF a THE SS TABE Lut Ye DaeMaer 

(2) THE CONVERGENCE OF B1,TE1,DT1 IS DETERMINED BY THE PARAMETER 

'IT', USUALLY A REASONABLE DEGREE OF ACCUACY IS OBTAINED BY 

SETTING ‘IT'=3 ,HOWEVER, *IT'=1 WILL SAVE TIME FOR SURVEY 

CALCULATIONS. 

(3) IN CASES WHERE IT iS EXPECTED THAT TI ECG fepe << eo 

CONVERGENCE IS SOMETIMES SPEEDED BY WRITING THE SECOND TEMPERATURE 

EQUATION IN TERMS OF T1({1I) RATHER THAN DVL) Fabry (Gays (VE elt ie CL rs) 

ITINC=0 

NLC=0 

O=P-1 

Oat aa ee 

RO (1) =0.1D-12 

RO (2) =W 

A= (R/W) ** (1.0/7 (P-1.0) ) 

pO 99 I=1,10 - 

A=(((R¥ (A-1.0)) +W) /W) ** (1.0/P) 

CONTINUE 

U0 (1) =0.0D0 

U0 (2) =0.0D0 

DO 1 I=3,9 

UO (I) =0.CD0 

RO (I) =RO (1-1) #A*(RO(1-1-R0 (I-2)) 

CONTINUE 

RO (P+2) =2*RO (P+1-RO (P) 

KB=1.6E-12 

H=13.6 

AV=7.0E16 
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PI=3.1416 
MI=1.6725E-24 
C=3.0E+10 
WL2=3.17E+28 
(e ASSUMES INITIALLY THAT TE=TI=1EV. 
BC (P+1) =Bz 
DOe2"Y=17P 
DR (I) =RO (1+1-R0 (1) 
VIS (1) =DR(I) *VI 
TEO(1)=1.0D0 
TPO(I)= 
TP2(L)= 
TI0(1)=1.0D0 
DTO(1) =0.0D0 
BO (I) = 
NO (I) =AV*PO 
DN (1) Beri (ROLLE EE Z. O- (RO (I) **2.0)) 
Q0(1I)= 
2 senah aaa 
T=0.0 
CALL REST (TIO,T,P,NTEST, NREC,W,R, DN, TPO) 
WRITE (6,101) PO,B2,T 
WRITE (6,102) 
WRITE (6,103) 
White (O77 104) (1,R0 (Tt) ,RO (1+1), UO(I), NO (1), TEO(ZT) ,DTO (1) ,20 (1), 
1I=1,P) 
101 FORMAT (1H1,' FILLING PRESSURE IS*',E10.3,10X," AXIAL MAGNETIC FIEL 
Tota pe leo! TIME SUS pits sos! SECS §/7/) : 
102 FORMAT ( ' THE INITIAL CONDITIONS ARE'/) 
AOS rOGNAT (( ! Tt OK, ROL) yk ye noite 1) ', 5X, oo COUT) "oro NO (te, 
Tike TEOCL) Y 6X," DIO (1). ",6X," BO (1) *7) 
104 FORMAT (14,7D13.4) 
J=0 
L=1 
M=C 
NLOOP=0 
NPRT1=NPRT 
ITEST=0 
CALL FLUX (P,BO,R0) 
4 CONTINUE 
NLOO P=NLOOP+ 1 
DIN=DT 
CICALCOLATE TIMESTEP *DT 
CAEL TAIN B(P DI; Liz,LTEST, KB, TIO, MI,O0, PP; PI) 
(e TEMP,DENSITY, AND MAGNETIC FIELD AT T+DT/2 INITIALIZED TO VALUES AT T. 
CALL ENERGY (DT,P,0) 
DO3 I=1,P 
RC (1) =(RO (1) +RO (I+1))/2.0 
TE2 (I) =TEO (I) 
DT2 (1) =DTO (I) 
TiAl) =L10 (1) 
N2 (I) =NO (I) 
B2 (1) =BO (I) 


3 CONTINUE 
Cc FIRST APPROX. FOR ELECTRON ION THERMALIZATION TIME + ELECTRICAL 
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CONDUCTIVITIES AT T+DT/2. (VALUES AT TIME T USED AS AN APPROX.) 


_ CALL ETCOND(KE2, P) 


CALL ITCOND(KI2,TI2,P) 

CALL EITT (TEQ2, 2) 

CALL ECOND(KR2,P) 

CALCULATION OF VARIABLES TE1,T11,B1,R1,01,N1,NE1, AT TIME T+OI: 
CALCULATION CF PONDEROMOTIVE PRESSURE 

CALL LHEAT(T,RC,HRATIO,P,1,0,T12,NCYCH,J,IT,DT) 

DOMUUSIT=1) 

PP (I) =3.18E+ C9*NOQ (I) *LIN (1) / (C¥WL2) 

CONTINUE 

SOLUTION FOR THE FLUID VELOCITIES AT T+.5DT: 

SOLUTION FOR THE SHELL BOUNDRIES AT T+DT: 

U1(1) =0.0D0 

U1 (Q) =U0 (Q) 

R1 (1) =RO (1) 

R1 (Q) =RO (Q) 

WILL DEFINE A NEW TIMESTEP "DTN* WHICH IS AN AVERAGE OF THE OLD 

TIMESTEP USED AND THE NEW TIMESTEP TO BE USED, OTN WILL BE USED 

FOR CALCULATION OF U1(I) ONLY. FOR FIRST TIMESTEP DTN=DT , 
DIN=(DTNtDT) /2.0 ; 

LPelT EST. EQ, 0) DIN=DT 

ITEST=1 

DO 5 I=2,P 

U1 (1) =U0 (I- (4. ODO *DTN* ( (KB*NO (L) * (TIO (LI) +TEO (1) ) #MI*QO (1) +PP (I)) 
1-(KB*¥NO (1-1) * (TIO (I-1) #TEO (1-1) ) *MI*QO0 (I-1) +PP(I- EBV AMES paar 
se 1) ) * (NO (1-1) +NO (I) ) *MZ) ) 

3- (DIN*2.0* (BO (I-1) +BO (L) )*(1.0/(4.0*PI)) * (BO (1-89 (I-1))/ 
4 ( (RC (I+1-RO (i-1)) * (NO (1-1) +NO(Z) ) *MI) ) 

R1 (I) =RO (1) + (DT* (01 (I) )) 

CONTINUE 

R1(2+2) =2*R1 (P+1-R1 (2) 

SOLUTION FOR THE NEW DENSITIES AT TIME T+tDT: 

CALCULATION “OF VISCOSITY AT TIAE® Te, SDT: 

DO 73 I=1,P 

N1 (I) =DN (I) / (( (81 (1*1) *R1 (1+1) -(R1 (1) *B1 (I) ) ) ¥PL) 

N2 (1)=(N1(1L)+NO(1)) 72.0 

CONTINUE 


CALCULATES DIFFERENCE QUANTITIES NEEDED FOR CALCULATING B1,TE1,DT1 


R4 (1) =RO (1) # KO (2) +R1(1) +R1(2) 

RC (1) =R4 (1) 74.0 

DR1(1) =RO (2) +R1(2- (RO (1) +R1(1)) 
DR1(P+1)=(KO (P+2) +R1 (P+Z) -(RO (Pt) +R1 (P+1)) 
DR2(1) =RO (3) +R1(3- (RO (1) +81 (1)) 

RC (P) =R 

QO (P) =0. 

DO 74 I=2,P 

DR1 (I) = (RO (+1) #R1 (1+1) - (RC (I) #R1 (Z)) 

DR2 (1) =(R1(1+2) +RO (1+2) - (R1 (1) + RO (1) ) 
DR3(1) = (R1(1+1) +RC (1+1) ~ (K1 (1-1) + RO (I-17) ) 
R4 (I) =RO (1+1) +R1(1+1) +RO (1) +R1 (1) 

PRNE NE Tey AOR ee OEE USED FOR NEXT TIME INTERVAL AS WELL. 


eynreet LT. U4 (1)) QO (I) =0.25¥*VIS (I) *VIS (I) *N2 (1) * ( (U1 (1+ 101 (I) 
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1) /DR1(I)) **2 

DX2 (1) =(R1(L+1) #RO (1+1)) *01 (1+#1) / (RY (1) *DR1 (I) ) 

DX3 (1) = (R1 (1) +RO (1)) *U1 (1) / (R4 (I) ¥DR1 (I) ) 

DX4 (I) = (R1(1+1) #RO (1+1)) /(R4 (1) *DR2 (L) *DR1 (1) ) 

DXS (1) = (R1 (I) +RO (1) )/(R4 (1) *DR3(1) *DR1(1)) 

DX6 (I) =(R1(L+1) +#RC (1+1)) / (RY (L) *DR2 (I) *DR1 (I) * N2(1)) 

DK7 (I) =(R1(I) + RO (I)) /(R4 (1) *DR1 (1) *DR3 (1) * N2(I)) 

BBB(I)= (1.0D0/DT) + ((2.0D0)* (DX2(1-DX3 (1))) 
DDB(1)=(BO(I) /DT) - ((2.0D0) *B0 (1) * (DX2 (1-DX3 (I))) 

BBI (1) =(1.0D0/DT) + ((4.0D0/3.0D0) * (DX2(I-DX3 (I) )) 

DDI(I)= (DTC (I)/DT) - ((4.0D0/3.0D0) *DTO (1) * (DX2(1-DX3 (1))) 
1+ ((8.0D0/3.0D0) *M1*Q9 (1) * (DX2 (I-DX3(I))/ N2(I)) 
BBE(1)=(1.0D0/DT) + ((4.0D0/3.0D0) * (DX2(I-DX3 (I))) 

DDE(I)= (TEO(I)/DT) - ((4.0D0/3.0D0) *TEO (1) * (DX2 (1-DX3(1))) 
CONTINUE 

T2=T+0.5¥*DT 

ITERATION LOOP FOR B1,TE1,T1I1 STARTS: 

J=-1 

CONTINUE 

J=JI+1 

CALL LHEAT(T2,R8C,HRATIO,P,2,1,T12,NCYCH, J, IT, DT) 

SOLUTION FOR THE MAGNETIC FIELD AT TIME T+tDT: 

EQUATION FOR THE MAGNETIC FIELD WILL BE WRITTEN IN THE FORM: 
-AA(I) B1(I+1) + BB(I)B1(1) - CC(I)B1(I-1) = DD (I) 

KR2 (P+1) =KR2 (P) 

DOn8 T=2;, 2 

RCL) = (DX4 (I) * (KR2(1+1)+KR2 (I) )/PT) 
BB(1)= BBB(L)+((1.0D0/PI) * (DX4 (I) * (KR2 (I+1) tKR2(1) ) +DX5 (I) * 
1(KR2 (1) +KR2(1I-1)))) 

Ce (tT) <a (DX5 (1) * (KR2 (1) *K R2 (1-1) ) /PI) 

DD (1) =DDB{1) +((1.0D0/P1) * (DX4 (1) * (KK2 (I+ 1) +KR2 (1) ) *( BO (I+1-BO0 (I) ) 
1-DX5 (L) * (KR2 (1) +KR2 (1-1) ) * (BO (1-BO (1-1)) )) 


8 CONTINUE 


9 


10 


CALCULATES MAGNETIC FIELD IN CENTRAL CORE: 
GG=2.0* (RO (2) +R1(2)) *U1 (2) / (DR1 (1) *B4 (1) ) #8.0* (RO (2) +R1 (2) ) * (KR2 ( 
12) #KR2(1))/(12.566*DK2 (1) ¥DR1(1) *R4(1))+1.0/DT 

EE=8,0*(RO(Z) +R1(2)) * (KR2(1) tKR2 (2)) / (DR2 (1) #12. 566*DR1 (1) *R4 (1) ) 
FF=8.0* (RO (2)+R1 (2)) * (KR2 (1) +KR2Z (2) ) ¥ (BO (2-BO (1) ) 7 (12. 566*DR1 (1) 
1*K4 (1) *DR2(1) -2.0*BO (1) * (RO (2) +R1 (2) ) *U1 (2) 7 (R4 (1) ¥DR1 (1) ) +B0 (1) / 
2DT 

E(1) =EE/GG 

F (1) =FF/6GG 

DO 9 I=2,P 

E (1) =AA (I) /( BB (I~ (CC (1) *E (I-1))) 

F (1) = (DD (I) + (CC (1) *F (1-1) )) / (BB (> (CC (I) *E (1-1) }) 

CONTINUE 

B1(P+1) =F (P) /(1-0-E(P)) 

B1 (P) =E(P) *B1(P+1) +F (P) 

DO 10 I=1,0 

B1 (P-L) =E (P-1) *B1 (Q-I) +F (2-1) 

CONTINUE 

Saini ol FOR ELECTRON TEMPERATURE AT TIME T+*DT; 

EQUATION FOR TE(T+DT) WILL BE WRITTEN IN THE FORM: 

~AA(I)TE1(I#1) + BB(1)TE1(1) - CC(1)TE1(1-1) = DD (1) 

pO 11 I=2,0 
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BE USES (EP Fara Bee ey ee aL) Pe) 
- . 73.0DO) * ((KE2 (I+1) + * 
1KE2 (I=1)) *DX7 (1) }) yA (I+4) #KE2 (I) ) *DX6 (1) + (KE2 (I) + 
CC (I) = (8.00073.0D0) * (KE2 (1) #KE2 (1-1) ) *DX7(I) 
AJUL (IL) =((2. OFKR2(1)) /(3.0*KB¥N2 (1) * (4.0 *PI) **2)) 
1*((((B1 (L+#1) +B0 (1+1) - (81 (1-1) #B9 (1-1) ) ) 7DR3 (I) ) **2) 
DD (I) =DDE(1) +( (8.0D0/3.0D0) * ( (KE2 (I+1) #KE2 (1) ) ¥DX6 (I) * (TEO (I#1- 
ITEO (1) ~ (KE2 (1) *KE2 (I-41) ) *DX7 (1) * (TEO (I-TEO (I-1)))) 
2- (DT2 (I) /TEQ2 (1)) #HEAT (I) tHJUL (I) 
11 CONTINUE 
CALCULATES TEMPERATURE RISE OF CENTRAL CORE: 
GG=1. 333% (RO (2) +R1 (2) ) #01 (2) / (DRI (1) *B4 (1) ) #5. 333% (RO (2) +R1 (2) ) * 
T(KE2 (2) #KE2(1) )/ ((N1 (1) #NO (1) ) *DR2 (1) ¥DRY (1) *R4(1)) #1.0702 
EE=5.333* (KO (2) +R1 (2) ) * (KEZ (1) #KE2 (2) ) /(DR2(1) *(N1 (1) #NO(1)) 
1*DR1 (1) *k4(1)) 
FF=5,333 * (RO (2) #R1 (2) ) * (KEZ (1) #K E2 (2) ) * (TEO (2-TEO(1)) 
FP=5.333* (RO (2) +R1 (2) ) * (KEZ(1) #KE2Z (2) ) ® (TEO (2-TEO (1)) 7 ((N1 (1) + NO 
1(1)) *DR1 (1) *R4 (1) *DR2 (1) -1,333*T EO (1) * (RO (2) +R1(2) ) *U1 (2) 7 (R4 (1) 
2¥*DR1(1)) + (TEO (1) /DT- (DT2 (1) /TEQ2 (1)) 
3+ (8.0/(3.0* ((4,0*PI1) **2.0) *KB)) * (KR2(1) #KR2 (2) ) * ( (BO (2) #B1 (2) -( 
4B0 (1) +B1(1)) )* ((BO (2) +B1 (2) - (BO (1) +B1 (1) )) 7 ((N1(1) #NO(1)) *DB2 (1) 
5*DR2(1)) +HEAT (1) 
E(1) =EE/GG 
F (1) =FF/GG 
Do #1292=2;10 
E (1) =AA (I) / (BB (I- (CC (I) *E(I-1))) 
F (I) = (DD(1) + (CC (1) *F (1-1) ) ) / (BB (I~ (CC (I) *B (I-1))) 
12 CONTINUE - 
TE1(P) =TE0 (Pe) 
DOMIs sl=1,0 
TE1(P-1) =(E(P2-1) *TE1 (Q-1)) +F (P-L) 
13 CONTINUE 
SOLUTION FOR RELATIVE TEMPERATURE, DT=TE-TI, AT TIME TtDT; 
DO 14 “I=2,0 
AA(I)= (8.0D0/3.0DO0) * (KI2 (I+1) #K12 (L) ) *DX6 (I) 
BB (I) =BBI (I) +((8.0D0/3.0D0) *({KI2 (I+1) #K 12 (1) ) *DX6 (1) + (K12 (I) + 
1KI2(1-1))*DX7(I))) +(1.0D0/TEQ2(I)) 
CC (L) =(8.9D0/3.0D0) * (KI2 (1) #Ki2 (1-1) ) *DX7 (I) 
DD (L) =DDI (1) + ((8.0D0/3.0D0) * ((KI2(1+1) #K 12 (1) ) *DX6 (I) *(DTO(It*1- 
4DT9 (I) - (KI2(1) #KI2 (1-1)) *DX7(1) * (DTO (I-DTO (I-1)))) 
2 -(1.0D0*DTO (1) /TEQ2 (1) ) +HEAT (1) +HJUL (1) 
7-((8.0D0/3.0D9) * ( (KI2 (1+1) #K12 (1) ) *DX6 (L) * ((TE1(1+1) +TEO (1+1)) 
B= (TE1 (I) +TEO (1) ) - (KI2 (1) #KI2 (1-1) ) *DX7 (1) * ( (TEO (1) +TE1 (1) - 
9 (TEO (I-1) +TE1(1I-1))))) 
A+((8.0D0/3.0D9) * ( (KE2 (I+1) *KE2 (1) ) *DX6 (i) * ((TE1(1+1) +PEO (I+1)) 
B- (TE1 (I) *TEO (I)) - (KE2 (I) *KE2 (1-1) ) *DX7 (L) * ( (TEO (1) +TE7 (1) - 
C (TEO (1-1) *TE1(I-1))))) 
14 CONTINUE 
CALCULATES RELATIVE TEMPERATURE IN CENTRAL CORE: 


Hi=oes 33 # (RO (2) +R1 (2) ) * (KI2 (1) #K12 (2) ) * ( (TEO (2) +TE1 (2) - 
1 (TEO (1) #TE1(1))) 7 ((N1 (1) #NO (1)) *DR1 (1) *R4 (1) *DR2 (1) ) 
Ti=5333 * (RO (2) #R1 (2) ) * (KE2 (1) +KE2 (2) ) * ( (PEO (2) #TE1 (2) - 


4 (TEO (41) #TE1(4)))/((N1 (1) #NC (1) ) *ORT (1) #84 (1) *DR2 (1)) 
5G=1635.3% (RO (2) #R1 (2) ) #U1 (2) / (DRI (1) #4 (A) ) +5. 333% (RO (2) #24 (2)) 
4* (KI2 (2) #K12 (1) )/ (0N1 (1) #NC (1) ) *DR2 (1) *DR1 (1) *R4(1)) +1.07DT+160/ 
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2TEQ2 (1) 

EE=5.333* (RO (2) +R1 (2) ) * (KI2 (1) *K12(2 DR2 * (N 

bean Me i (1) + (2))7( (1) * (N17 (1) #NO(1)) *DR1 
PF=5.333 * (RO (2) #R1 (2) ) * (KI2 (1) +KI2 (2) ) * (DTD (2-DTO (1)) 

TZ ((N1 (1) + NO (1) ) *DR1 (1) *RY (1) *DR2 (1) -1. 33. 3*DTO (1) * (RO (2) #81 (2)) * 

2U1 (2) /(R4 (1) *DR1(1)) #DTO (1) /DT-DTO (1) /TEQ2 (1) 

ts (8.0/7 (3.0* ((4,0*PI) **2.0) *KB)) *(KR2(1) +KR2(2)) * ( (BO (2) +B 1 (2) - ( 
SBC (1) +B1(1)) ) * ( (BC (2) #B1 (2) - (BO (1) +B1 (1) ) ) 7 ( (N17 (1) #NO (1) ) *DR2 (1) 


15 


16 


17 


20 


222 


6*DR2(1)) *HEAT(1-HH+II 

E (1) =EE/GG 

FP (1) =FF/GG 

DOMIa.5=2)50 

E(1) =AA (I) /(BB(I- (CC (I) *E (I-1))) 

F (1) =(DD(Z) + (CC(I) *F (I-1))) / (BB(I- (CC (I) *E(I-1))) 
CONTINUE 

DT1(P) =DTO (P) 

DT1(0) =E(0) *DT1(P) +F (0) 

DOIG T1=2 0 
DT1(P-1) =E(P-L) *DT1(Q-1) +F (P-I) 
CONTINUE 

CALCULATION OF TRANSPORT COEF AT TIME TtHDT/2. 
DOM a1, 2 

TE2 (1) =(TEO (I) +TE1(I)) /2.0 

DT2 (1) = (DTO (I) +DT1(1)) /2.0 

TI2(1) =TE2(1I-DT2 (1) 

B2 (1) = (BO (1) +81(1I))/2.0 

CONTINUE 

CALL ECOND (Kk&2,P) 

CALL ETCOND(KE2, P) 

CALGEETCOND UKI2,T12, 2) 

CALUSELTT (TEO2 /P) 

RETURNS TO CALCOLATE TE1,B1,DT1 ONCE MORE 
GENERALLY FOUR ITERATIONS WILL SUFFICE 
LE-s(J-1T) 77;20,20 

CONTINUE 

NOW HAS FINAL VALUES FOR B1,TE1, DT1 

IT INC=ITINC#+1 

T=T+DT 

RELABEL VARIABLES VARIABLE1 BECOMES VARIABLEO. 
RO (Q) =R7 (Q) 

RO (P+2) =R1(P +2) 

BO (P+1) =B1(P+1) 

UC (Q) =U1 (Q) 

DO 21 I=1,P 

RO (I) =R1 (1) 

UO (I) =U1 (1) 

NO (I) =N1(1) 

TEO(L)STE1 (L) 

DTO (1) =DT1 (1) 

TTO(Ly=CH4I(I-DT1I (1) 

BO (I) =B1(1) 

CONTINUE 

IF (L-PRTFR) 223,222,222 

CONTINUE 

CALL ENERGY (DT,P,1) 
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DOES O5e1=1,2 

POWER (I) =PABS (I) /(3.14¥*(R1(I+1) **2-R1 (I) **2)) 
505 CONTINUE . 

NERTI=NPRT1+1 

IF (NPRT1.LT.NPRT) GO TO 881 

NPRT1=0 

WRITE (6,105) 1, NLOOP 

105 FORMAT ($19, TINE S=0" 612.6," SEC", 
1° THIS IS THE *,I14,* TH TIMESTEP? /) 

WRITE (6,106) 

106 FORMAT(* *,* (1)*,1T10,"RO (I) ',T20,"UG (I) *,T31,*TEO (I) *,T42,"DI0 (I) 
1, T54,"NC (1) *,T65,"BO (1) ',T75,"PABS(1) ',T87,*PTH(1)*,197, "PKE(I) * 
2¢7108, *PB(I) *,1118, "HRATIO(I) '/) 

DO225 1=1,2 
WRITE (6,107) I, RO (I) ,UO(L) , TEO (1) ,DTO (I) , NO (I) ,BO (I), PABS (I) , PTH(I) 
1,PKE(I) ,PB(L) , HRATIO(1) 

225 CONTINUE 

107 FORMAT(" 272476 (1X,D10~3) ,5(1%,£1023)) 

WRITE (6,108) PABST,PTHT, PKET, BT 

108° FORMAT (#08) PABST® =!7810 53," = PTHTs=",E10.3," PKET =0,610.3 
Jel PP BTA=2,61033) 

CALCULATE AND PRINT OUT VALUES USED AS AN ENERGY CHECK ON PARALLEL 
TO PERPENDICULAR ENERGY TRANSFER. 
WRITE (6,501) 
501 FORMAT(*1"',* (I) PAR. TEMP POWER IN POWER OUT POWER ABS, TOTAL 
1 POWER E3"') - 
TW3M=0.0 
TW3M1=0.0 
TDPAR=0.0 
po’ 508 I=1,P 
W3MW=W3M (Z) *1.0E-07*DN (I) /N2 (I) 
W3MTW=W3M1 (1) *4.0E-07*DN (1) /N2 (I) 
DEPAR=0.8E-19*TPAKD (1) *DN (L) /DT 
TW3M=TW3M+W3MW 
TW3M1=TW3N1+ W3M1W 
TDPAR=TDPAR+DEPAR 
WRITE (6,502) I,TPF (1) ,W3MW, W3M1W, DEPAR, POWER(I) , £322 (1) 
508 CONTINUE . 
502 FORMAT(" ',14,6E11.3) 
WRITE (6,507) 

507 FORMAT('O',* TOTAL VALUES FOR POWER ABSORBED AND RELEASED BY PARALL 
1EL TEMP. °) 

WRITE (6,506) TW3M,TW3M1,TDPAR 

506 FORMAT(* *,* POWER IN=",E10.3,'° POWER OUT=*,E10.3, 
an POWER ABSORBED=',E10. 3) 

881 CONTINUE 

L=0 
DATA READ ONTO A TAPE IF ITAPE= 1. 
IF (ITAPE.EQ.0) GO TO 223 
| CALL TPDATA 
223 CONTINUE 
IF (T-TF) 23,23,24 
23 L=L+t1 
GO TO 4 
24 CONTINUE 
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IF (ITAPE.EQ.1) ENDFILE3 

CALL FLUX(P, BO,RO) 

RETURN 

END 

SUBROUTINE LHEAT(T,RC,HRATIO, NCELL, NLASER, LHRT ,T12,NCYCH,J,1T, Df) 
REAL LIN(100),HKEAT(10C),RC (100) , TE2 (100) ,N2 (100) , 82 (100) , HRATIO (10 
19) ,LMDA,K3,PWR1 (100) ,P¥R2(100) ,w# 3H (100), TI2 (100) 
1, TPF (109) ,TPO (100) ,TP2 (100) ,W3M1 (100) , E322 (100) , TPARD (100) 
CONMON/COMTHZ/TE2,N2,62 /COMLAS/LIN,HEAT/COMUDIF/WAVE1,WAVE2, PWR1 
1,PWR2 /PARAL/TPO,TPF,TP2,E322,W3M1,TPARD, W3M 

LHEAT CALCULATES LASER INTENSITY LIN(I) (ERGS/CM2.SEC) (KEQUIKED 
FOR PONDERMOTIVE PRESSURE AND HEATING RATE/PARTICLE (uV/SEC) 
IHRT=1 IF HEATING RATES ARE TO BE CALCULATED. 


NCELL = NO. OF CELLS USED IN DIFFERENCE SHCHEME {=P IN NEAL'S NOT.) 


NLASER = NO. OF LASERS USED IN HEATING, 

HRATIO(I) = DIFF FREQ, HEAT/INVERSE BREMSSTRAHLUNG HEATING BY 
INCIDENT BEAMS, 
DEFN OF LASER PARAMETERS, 

WAVE1, WAVE2 - WAVELENGTH (MICRONS) $1, S2 - SPOT SIZE(CM), 
PD = PULSE DURATION (SEC) 

WAVE1=9.6 

WAVE2=10.6 

S$1=0.045 

S2=0.045 

P1=1.25E+9 

P2=1.25E+9 

PD=40.,0E-09 

LASER AMPLITUDE VS. TIME (AF= AMPLITUDE FACTOR) 
AF=2.0-2.0*T/PD 

TF (T.LE.PD/2.) AF=2.0*T/PD 

P1A=P1*AF 

P2A=P2*AF 

A1=2.0*S1*S1 

A2=2.0*S2*S2 

SPATIAL DISTRIBUTION OF LASER POWER. 

PWR1(I) = INTENSITY (W/CM2) 

LIN = TOTAL LASER (1+Z) INTENSITY (ERGS/CM2.SEC) 

DO2 I=1,NCELL 

PWR1 (I) =P1A/ (3.142*A1) *EXP (-RC (I) *RC (I) /A1) 

LIN (1) =1.0E+07*PWR1 (I) 

GO TO (2,1) ,NLASER 

PWR2 (I) =P2A/ (3.142¥*AZ) *EXP (-RC (1) *RC (I) /A2) 

LIN (1) =LIN (1) +1. E+07*P WR 2 (I) 

CONTINUE 

IF (IHRT.NE.1)GO TO 3 

CALCULATION OF HEATING RATE FROM INCIDENT LASER BEAMS VIA INVERSE 
BREMSSTRAHLUNG ABSORPTION. HEAT IS IN ERGS/CM3.SEC. 

Z=1.0 

DOS I=1,NCELL 

TX=TE2 (I) **1.5 

Ws One tees bs Loe 

oh (TE? (L)'.uts27. LMDA=2. 

HEAT (x) 20.67 £231" ((NZ (1) #WAVE1) #¥2) #Z*ALOG (LMDA) ¥PWR1 (I) /TX 
GO TO (5,4) ,NLASER 


HEAT (I) =HEAT (I) +8. 67E-31/TX* ((N2 (I) *WAVE2) **2) *2*A LOG (LMDA) *PWR2 (I 
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W3M1 (I) =0.0 

CONTINUE 

IF(NLASER,NE.2)GO TO 7 

IF(NCYCH.EQ.0) GO TO 10 

CALL DFHH (NCELL) 

DOS V1l=1, NCELL 

TCOL= (1.8E+0 4 (TE2 (I) +0. 333333*TP2 (I) ) **1.5) /N2 (I) 

rely etieo ise( O-ocsebt/zc 

P(I)= (1.0-0.5*DT/TCOL) +2.0*W3! * »6E-12* 

Wie oe oe sean TeoLt ua ) O*W3M (I) ¥DT/(1.6E-12*N2 (L))) 
TP2 (I) =0.5* (TPO (1) +TPF(I)) 

W3M1 (1) =1.60E-12*TP2 (1) *¥N2 (1) /(2.0*TCOL) 

GO TO 11 

TPF(1)=0.0 

TP2(1L) =0.0 

W341 (1) =W3M (I) +0.5*1, 6£-12*T PO (L) *N2 (I) 

CONTINUE . 

IF(J-IT) 12,13,13 

CONTINUE 

DOes T=17NCELL 

TPARD (1) =TPF (I-TPO (TI) 

TPO(L) =TPF (I) 

CONTINUE 

CONTINUE 

DO6 I=1,NCELL 

IF (HEAT (L).LT.1,0E-65) HEAT (I) =1.0E-65 

HRATIO (I) =W3M(1) /HEAT (1) 

HEAT (L) =HEAT (1) +W3M1 (I) 

CONTINUE 

2/3. HEATING RATE (EV.SEC)/PARTICLE VALUE STORED IN HEAT(I1). 
DOS I=1,NCELL 

HEAT (1) =0.6667*0,.625E+12 *HEAT (I) /N2(1) 

CONTINUE 

RETURN 

END 

SUBROUTINE DFHH(NCELL) 

DIFF. FREQ. HARMONIC HEATING VERSION 2 

ROUTINE CALCULATES HEATING RATE DUE TO MIXING OF TWO ANTIPARALLEL 
LASERS. 

COMPLE XP ESZ, FO ,U7/ (0.0, 1.0)7 

REAL T(100),N(10C) ,B(10C) ,PWR1 (100) ,PWR2 (100) , ¥34(100) ,K3,KZ,K 
1,FOSPL(12,5) ,TP22(100) ,E32Z2 (100) , TPO (100) , TPF (100) ,W3M1 (100), 
1TPARD (100) 

COMMON /COMTH2/T,N, B/COMDIF/WAVE1, WAVE2,PWR1,PWR2 
1/PARAL/TPEO,TEF,TP2, £322, W5N1,T PARD, W3M 

DATA ePOSPL/O e210 5,0c8, 100, 1ele lady lsle2sO yee de 300 sand te 0, 
1 .3895,0.7199,1.064,1.076,1.015,0.913,9. 7451,0.6027,0.4462,0.3565 
2,0.2992, 0.2587, 1-903 412299,0.5474,-0.2165,-0.4162, -0.5635,-0.532, 
3-0.4114,-0.2316,-0613934-0.0931,0.0,0.9,0.987,-344911,-0.3284, 
4-0.6701,-2-0 6606, 021709, 0.231,0.1287,0.5593E-1,0403642,0.0,1.997 
59 -4.976 450272 ,-04 569641 007,04 2633,0.06673,-0.06821,-0.0485, 
6-0.01301,-0.02428,0.0/ 

CALCULATION OF BEAT FREQ. LASER HEATING, 

PARAMETERS FOR THE DIFF. FREQ. WAVE (ANTIPARALLEL BEAMS) 
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K3=6.283E+04* (1,0 /WAVE1+1.0/WAVE2) 
W3=7.885E+15* (1.0/WAVE1-1.0/WAVE2) 
KZ=ABS (K3) 

DOT T1=4,NCELL 
WPE=5.63E+04*SORT(N(II)) 
VOE=5.93E+07* SQRT ((T (IL) +TP2(II))) 
VET=2,53E-O5*N (II) /((T(I1) #0.333333*TP2 (IL) ) **1.5) 
K=KZ*VOE/ (1.4142 *WPE) 

AQ=W3/ (KZ*VOE) 

AQE=A0 

IF (ACE.GT.10.) AOE=10.0 

AO ES=AQO E*AQE 


SPLINE APPROXIMATION FOR FO 

IF (A0.GT.3.25)G0 TO 12 

bO7 T=1, 12 

IF (AC.LT.FOSPL(I,1))GO TO 8 
7 CONTINUE 
8 I=I-1 

DA=A0Q-FOSPL(I,1) 

FO=FOSPL(1I,2) +FOSPL (1,3) *DA+FOSPL(I,4) *DA*DA+FOSPL (1,5) *DA¥*¥*3 

FC=J *FO . 

GO T0814 . 
ee FO=J/A0*(1.04+0.5/(AO*AG) +0. 75/7 (A0¥*4) +1.875/(A0¥*6) +6.565/ (AO**8) ) 
11 DAMP=AMAX1 (VEI/(2.0*W3*A0**3) , 1. 772*EXP (-AQES) } 

FO=FO +#DAMP 

E3Z=3*3 .0E-16¥*KZ*WAVE1*WAVE2* (1,0 +3 *AO*FO) / (K¥K41.0+0*AO*FO | 

1-2. 72E-04/ (AO*A0) ) *1.0343 
E322 (11) = (CABS (E32) ) *SQRT(PWR1 (IL) *PWR2 (II) ) 2 
W3M (IL) =WPE*WPE*AO*#A0/ (3. 142*KZ*VOE) *REAL (FO) * (CABS (E3Z) **2) 
1*PWR1 (II) *PWR2 (II) 

1 CONTINUE 

RETUBN 

END 

SUBROUTINE REST(TIC,T,P2,1TEST,NREC,W2,R,DN,TLO) 

THIS ROUTINE WILL READ DATA FROM SEQ. PILE AND PASS IT BACK TO 

ALAS2. 

AS INITAL CONDITIONS. 

IF ITEST=0 NOTHING IS DONE, 

IF ITEST=1 ROUTINE RESTARTS WITH INITIAL CONDITIONS AS ON FILE. 

IF ITEST=2 ROUTINE REORGANIZES DATA TO CHANGE NUMBER OF SHELLS, 

REAL*8 RO (101) ,U0(101) ,TEO (100) ,DTG (100) , TIO (100) , 80 (100), 

1N0(10C) ,DBLE,DN(100) 

REAL RS (102) , TES (102) , TIS (102) ,NS (102) ,¥S (102) ,BS(102), 

1RCS (100), RCN (100) , TPS(102) , TPO (100) 

INTEGER P1,P2,Q0 

COMMON/CONTHO/TEC, DTO, NO,BO/COMHYO/U0, RO 

TPUITESTIEO; 0) GO TO" 1 

DO 2 I=1,NREC 

READ (4) T, TF, P0,R,BZ,P1,RS,TES,TIS,NS,VS, BS, TPS 

2 CONTINUE 

TEUITEST. EQ. 2) GO TO 9 

pows i=17P1 

RC (I) =DBLE(RS(1)) 

UO (I) =DBLE(VS(1)) 
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TEO (I) =DBLE (TES (I) ) 

T10(1) =DBLE(TIS(I)) 

TPO (1) =TPS(I) 

BO (I) =DBLE(BS(1I)) 

NO (1) =DBLE(NS (I) ) 

DTO (I) =TEO(I-T10 (1) 

CONTINUE 

BO (P1+1) =DBLE(BS (P1)) 

RO (P1+1) =DLE(RS (P1+1) ) 

RO (P1+2) =DBLE(2.0*RS(P14+1-RS (P1) ) 
UO (P1+1) =DBLE(VS (P1+#1) ) 

DO 14 I=1,P1 

DN (I) =NO (I) *3.1416* ((RO (i+1) **2- (RC(L) **2) ) 
CONTINUE 

GO TO 1 

RQ (1) =0.1D-12 

RO (2) =w2 

A= (R/W2) ** (1.0/(P2-1.0)) 

DO 4Y I=1,10 

A=(((R¥ (A-1.0) )+W2) /W2) ** (1.07/22) 
CONTINUE 

Q=P2+1 

DOP Seel=37,.0 

RO (I) =RO (1-1) +A* (RO (I-1-RO (I-2)) 
CONTINUE 

RO (Q+1) =2.0* KO (P2+1-R0 (P2) 

DO 6 I=1,P1 

RCS (1) = (RS (1) +RS (I+1)) /2.0 

DO 7. 1=1,e2 

RCN (L) = (RO (1) +RO (I+1)) /2.0 

N1=1 


DO 10 J=N,P1 
N1=J 
IF (RCN(I).LT.RCS{J)) GO TO 11 
CONTINUE 
ONLY COMES HERE IF RCN(I) .GT. RCS(P1) 
UC (I) =DBLE(VS(P1)) 
TPO(1) =TPS(P1) 
TEO(1) =DBLE(TES (P1)) 
T1O(1) =DBLE(TIS (P1)) 
BQ (I) =DBLE(BS(P1)) 
NO (I) =DBLE(NS(21)) 
DTO(1) =TEO (I-TIO (1) 
GO TO 8 
IF(N1.EQ.1) GO TO We Pen ety 
I) IS BETWEEN RCS (N1- 
BO (4) =BS (1-1) #( (RCW (I-RCS (N 1-1) ) 7 (RCS (NT-RCS (N1~1)) ) 
* (BS (N1-BS (N1-1 
TF h tr pans 11-1) + (RCM (T“RCS (N1-1) ) 7 (HES (H4-RCS(N1-4))) 
x i1-V5(N1-1 
ree Uy eans N4=4) © ( (RCN (Z-RCS (N1-1) )/ (RCS (N1-RCS(N1-1))) 
% N1-TES (N1-1 
weed ea hetas (1-1) © CURCN (Z-RCS (N1-1)) / (BCS (N1-RCS (1-1) )) 
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1* (TIS (N1-TIS (N1-1)) 
TPO(L) =TPS(N1-1) + ( (RCN (I-RCS (N1- i1- i1- 
18 (22S (W185 (1-1) ( ( 1))7 (8CS (N1-RCS (N1-1))) 
NO (1) =NS (N1-1) #( (RCN (1-RCS(N1-1 cs - - 
1# (NS (N1-NS (N4-4)) ( ( )) 7 (RCS (N1-RCS (N1-1)) ) 
DTO(1) =TEO(1-T10(I) 
GO TO 8 
12 U0 (I) =DBLE(VS(1)) 
TEQ(I) =DBLE (TES (1) ) 
TPO(1) =TPS (I) 
T1O(1) =DBLE(TIS(1)) 
BO (I) =DBL&(BS (1) ) 
NO (1) =DBLE(NS(1)) 
DTO (1) =TEO(I-T10 (1) 
8 CONTINUE 
BO (P2+1) =B0 (P2) 
UO (P2+1) =0.0D0 
DO 13 1=1,P2 
DN (I) =NO (I) *3.1416*( (RO (1+1) **2.0-(RO (I) **2.0) ) 
13 CONTINUE 
1 RETURN 
END 
SUBROUTINE TIME (PP), TSF, TI2 LTEST KS, TI10,N2, 00, PP, P1L) 
REAL*8 TE1(100) ,07T1(100) ,N1 (100) ,5b1 (100) ,01(161) ,R1(101), 
180 (101) ,U9 (101) ,TEO (100) , DTO (100) ,B0 (100) , NO (7100) 
2,TIO (100) ,g0 (106) , PP (109) , DABS 
REAL TE2(160),N2(100),T1I2(100) ,82 (109) ,KB,NI,AA (101) ,T2 (6) 
INTEGER P,Q 
COMMON/CONTH1I/TE1,DT1,N1,B1/COMHY1/01,R1/CONTH2/TE2,N2 
1,B2/COMHYO/UC,RO/COMTHO/TEO ,DTO, NO, BO 
SUBROUTINE CALCULATES TIME INCREMENTS FOR NEXY STEP IN MHD CODE 
P=NO. OF CYLINDRICAL SHELLS USED 
FIRST CALCULATE TIME STEP FOR EQNS. TO BE STABLE. 
TSF=FINAL TIME STEP TO BE. USED 
TSI=TIME STEP USED IN PREVIOUS LOOP 
CS2=(SPEED OF SOUND) **2 
CA2=(ALFVEN VELOCITY) **2 
TZ (1) =1000.0 
TZ (2) =1009.0 
TZ (3) =1000.0 
TZ (4) =1000.0 
TZ (5) =1000.0 
TZ (6) =1000.0 
A0=0.8 
A1=0.2 
A2=1.0 
A3=1.0 
A4¥=10 
A5=20.0 
A6=0.2 
TSI=TSF 
TSF=1.0 
IF (ITEST.EQ.0) GO TO 601 
GO TO 602 
601 CONTINUE 
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DOnS eis 1, P 
TE1 (1) =TEO (I) 
DT1(1) =DTO(I) 
N1 (I) =NO (I) 
B1 (I) =B0 (I) 
R1 (L) =RO (I) 
3 CONTINUE 
AQ=0.1 
R1(P+1) =RO (P +1) 
602 CONTINUE 
ce LINIT TS SO STABILITY CONDITION IS SATISFIED, 
DO. 1, L=1,P 
PRESS=N1(I) * (1, 380-16) * (TE1 (I) + (TE1 (I-DT1 (I))) 
RHO=(1.670-24) *N1 (I) 
CS2= (5.0/3.0) *PRESS/RHO 
CA2= (B1 (I) **2) 7(4.0*3.14159*KHO) 
C=SQRT(CS2+CA2) 
TS=A0*(R1(1+1-R1(1)) /C 
TPCT ol lial oi 1) eco iets 
1 CONTINUE 
EP O(LTEST,EO.0) GO. T0_99 
Q=P+1 
C CALCULATE AA(I) =ACCELERATION OF THE I*TH SHELL BOUNDARY, 
AA(1) =0.0 
AA (P+1) =0.0 
DO 6 I=2,11 
AA (L) =- (4.0D0* ( (KB¥NO (I) * (TIO (I) +TEO (1) ) +*MI*QO (1) +PP (1) ) 
1-(KB¥*NO (1-1) * (TIO (1-1) +TEO (1-1) ) #M1*Q0 (1-1) +PP (I-1)) ) /( (RO (1#1- 
2RO (1-1) ) * (NO (1-1) #NO (I) ) MZ) ) 
3-(2.0* (BO (1-1) +BO (1) ) * (1.0/7 (4.0*PTI) ) * (BO (1-BO (1-1) )/ 
4 ( (RC. (1+1-RO (1-1) ) * (NY (1-1) +NO (1) ) *MI) ) 
IF (ABS (AA(I)).LT. 1.0E-60) AA(I) =1.0E-60 
6 CONTINUE 
DO 2 I=1,10 
C LIMIT TS SO SHELL THICKNESS IS NOT INCREASED BY MORE THAN A 
CEPACTORSOF et A1*, 
RI=A1* (R1 (1+ 1-81 (I)) 
U=01(1+101 (1) + (TSI* (AA (I+1-AA(I) ) /2.0) 
AB=(AA(I+1-AA(I))/2.0 
AU=ABS (U) 
AAB=ABS (AB) 
IF (ABS (U) .LT.1.90E-30) U=1.0E-30 
LF (AAB.LT.1.0E-50) AAB=1. 0E-50 
IF (((4.O*AAB*RI) /(U**Z)) -GT.1.0E-02) GO TO 7 
TS=RI/AU 
GO TO 8 
7 DDIS=AU**2-4,0*AAB*RI 
SDIS=AU**2+4 .O*AAB¥RI 
TF (((J+GT.0.0) «AND. (AB.GT.0,0) ) «ORs ((0.LT.0.0) .AND. (AB. LT.0.0))) 
=(-AUtSQKT (SDIS) )/ (2-9 *AAB) 
Te C1 (U.G2. 0.6) AND; (AdcL?.0.0) -ANDs (DOIS-67.0.0)) 08. ( (U.LT.0.0) 
1. AND. (AB.GT.0.0) .AND. (DDIS.GT.0.0))) TS= (AU-SQRT (DDIS) )/(2.0*AAB) 
LF (((U-+GT.0.C) -AND. (AB.LT-0.0) sAND. (DDIS.LT.0.0)). OR, ( (U.LT.0.0) 
1. AND. (AB.GT.0.0) «AND. (DDIS.LT.0.0))) TS=(AUFSQRT(SDIS) )/(2.0*AAB) 


8 IF(TS.LT.T2(2)) T2(2)=TS 


85 


cat) hd ees 
=f r2POR) PLB t- Gaiters 
, an 
Ni ttosI ere 
wet | 
ave 7 ; 
A WANS MOK TE AeSsse Il: Bes ES es 


or 
(WwS\.4 (issih e tp iy 
&.2y 


et “ ial 
5 OT O68 (40-30, Pog ; 


grist itn FOR GaSe TIM AD, ( 10. 0475 « 


C LIMIT TS SO CELL ELECTRON AND ION TEMP, DOES NOT CHANCE BY HORE 
C THAN A FACTOR OF *A2* AND 'A3* RESPECTIVELY. 
T=2.0* (ABS (SNGL (TE1 (1) -TE2(I))) 
reel Ut. OE =50)0 T=1,. 108450 
TS=(A2* TET (1) *TS1) sr 
ERU(TSOLT . TZ:( 3) ) TZ (3) =TS 
MII =STEW(L-DTA (ry 
EES 2iO* (NBSP EI= 11 2°(1) )) 
IF (TL.LT.1.086-50) TI=1,0E-50 
TS=(A3*T11*TSI) /TI 
PRCTS UT S224) )11TZ (4) =TS 
2 CONTINUE 
C LIMIT TS SO VELOCITY DOES NOT CHANGE BY MORE THAN A FACTOR OF AS. 
DO 4) T=2,10 
IF (DABS (U1(I))..LT.1.0D-50) U1(I)=1.0D-50 
IF (DABS (U1 (I)) .LT.0.1D+06) A5=50.0 
IF (DABS (U1(I)) .LT.C.1D+04) AS=0.1D+#06/DABS (U1 (I)) 
AV EL=DABS (AA (I) /U1-(I) ) 
IF (AVEL.LT.1.0E-5C) AVEL=1.0E-50 
TS=A5/AVEL 
re( TS Vet25)) £2.(5)=TS 
4 CONTINUE 
C LIMIT TS SO RADIUS DOES NOT CHANGE BY MORE THAN A FACTOR OF A6,y 
DOS 2/410 
RI=A6*R1 (1) 
U=U1 (I) +(TSI*AA(I) /2.0) 
AB=AA (I) /2.0 
AU=ABS (U) 
AAB=ABS (AB) 
IF (ABS (U).LT.1.9E-30) U=1.0E-30 
IF (AAB.LT.1. 0E-59) AAB=1.0E-50 
LE (((4.0*AAB*R1) / (U**2)) .GT.1.0E-02) GO TO 17 
TS=Riyau 
GO TO 18 
17 DDIS=AU**2-4.0*AAB*RI 
SDIS=AU**2+4,0 *AAB*RI 
IF (((U.GT.0.0) AND. (A8.GT.0.0)).OR. ((U.LT.0.0) -AND. (AB.LT.0.0))) 
1TS=(-AU+SQRT (SDIS) )/(2.0*AAB) 
IF (((U.GT.0.0) AND. (AB. LT.0.0) .AND. (DDIS.GT.0.0)) 08. ((U.LT.0.0) 
1. AND. (AB.G7.0.0) -AND. (DDIS.GT.0.0))) TS=(AU-SQHT(DDIS) )/(2.0*AAB) 
IF (((U.GT.0.C) .AND. (AB.LT.0.9) -AND. (DDIS.LT.C.0)). OR. ( (U.LT.0.0) 
1.AND. (AB. GT.0.0) -AND. (DDIS.LT.0.0))) TS=(AU+SQRT (SDIS))/(2.0*AAB) 
18 LP (TSVLIJTZ(6)) [2 (6) =TS 
5 CONTINUE 
C ONLY ALLOW TS TO INCREASE BY FACTOR OF AY IN A SINGLE STEP 
TS=A4¥*TSI 
Pri TS AL heTSF), @SF=TS 
99 CONTINUE 
c SET TSF EQUAL TO SMALLIST OF T2(I). 
nowe20) 1-106 
IF (T2(1).LT. TSF) TSP=12(1) 
0 CONTINUE 
; WRITE(7,102) T2(1) ,T2(2) -TZ(3) -T2(4) ,T2(5) - TZ (6) 
WRITE(7,101) TSF 
102 FORMAT (*°0',6E15.3) 
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RETURN 

END 

SUBROUTINE ETCOND(KE,NCELL) 

SUBROUTINE CALCULATES ELECTRON THERMAL CONDUCTIVITY AT SHELL 
CENTERS. 

REAL = TE (100) ,N(10C) ,B (100) ,KE (100) 
COMMON/COMTH2/TE,N,B 

DO1 I=1,NCELL 

GA=2. 9E-O5#N (I) /TE (I) **1,5 

WCE=1,. Peet bt Se) 

KE (1) =2.556E+15*N(1I) *TE (I) / (GA¥(1.0+0.555* ( (WCE/GA) **2) ) ) 
CONTINUE 

RETURN 

END 

SUBROUTINE ITCOND(KI,TI, NCELL) 

SUBROOTINE CALCULATES THE ION THERMAL CONDUCTIVITY AT THE CELL 
CENTERS, 

REAL TE (100) , TI (100) ,N (100), 8 (100) ,K1(100) 
COMMON/COMTH2/TE,N,B 

DOA T= NCELL 

BE=1.414+0.175*(TI(I) /TE (I) ) **1.5 
GA=6.77E-07*N(I) /TI (I) **1.5 

WCI=9.56£E+03 *B (I) 

KE(E) 2S 98E+ 12*N (1) *TI (I) /(BE*GA* (1. 0+6,25* (WCL/ (BE*GA) ) **2) ) 
CONTINUE 

RETURN 

END: 

SUBROUTINE EITT(TEQ,NCELL) 

SUBROUTINE EITT CALCULATES ILON- ELECTRON THERMALIZATION TIME. 
REAL TE(100).,N (100) ,8(100) ,TEQ (100) ,LMDA 
COMMON/COMTH2/TE,N,B 

DOwd L=1,NCELL 

TH (TE (1) oL74736.0), GO TO:2 

LMDA=9.3170E+190*TE (I) /N (1) **0.5 

GO TO 3 

LMDA=1. 5554+10* (TE (I) **1.5) /N (I) **0.5 

TEQ(I) =3.16E+O8*TE (1) **1.5/ (N (1) ¥ALOG (LM DA) ) 

CONTINUE 

RETURN 

END 

SUBROUTINE ECOND(KR,NCELL) 

SUBROUTINE CALCULATES C**2 X RESISTIVITY (STATOHMS.CM3/SEC2) 
REAL TE(100),N(1C0),B (100) ,KK (100) ,LMDA 
COMMON/COMTH2/TE,N,B 

DOsmi=1,NCELL 

IF (TE (1)+LT.36.0)G0O TO 1 

LMDA=9.3160E+10*TE(I)/N (1) **0.5 

GO TO 2 

LMDA=1. 55E+10* (TE(1) **1. 5) /N (1) **0.5 

KR (1) =1.03E+O7*ALOG (LMDA) /TE (I) **1.5 

CONTINUE 

RETURN 

END 

SUBROUTINE ENERGY (DT, NCELL, LE) 
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ENERGY SERVES AS A CHECK ON THE SIMULATION PROBLEM. IN PARTICULAR 
THE RATE AT WHICH THERMAL, KINETIC ANC MAGNETIC ENERGY IS CHANGING 
IN A CELL IS CALCULATED AND CAN BE COMPARED TO THE RATE AT WHICH 
LASER ENERGY IS DEPOSITED IN THE PLASMA, 

REAL*8& UD (101) ,RO (101) , 726 (100) ,DT0(100) ,NO (100) ,B0 (100) ,DN(100) 
REAL ETHO (10C) ,ETH1 (100) ,EKEQ (100) , EKE1 (100) , EBO (100) ,EB1 (190) 
»PABS (1090) ,P5 (100) ,PTH (100) , PKE (100) ,HEAT (100) , LIN (100) 

COMMON /COMHYO/U0,RC/COMTHO/TEO,DTO,NO,BO/COMLAS/LIN, HEAT/COMEN/DN 
eETHC,EKEO,EB0,PABS,Pb, PTH, PKE, PABST, PTHT, PKET,PBT 

CALCULATION OF THERMAL,KINETIC AND MAGNETIC ENERGY IN EACH CELL, 
SET IE=0 IF THIS IS THE ONLY CALCULATION TO BE MADE. 
IF(IE.NE.0)GO TO 5 

DO6 I=1,NCELL 

ETHO (I) =2.4E-12*DN (1) * (2.0*TEO (L-DTO(L)) 

EKEO (I) =2.091E-25*DN(L) * (UO (1) +U0 (I+1) ) **2 

EBO (I) =BO (1) *BO (I) /8.0* (RO (L+1) *RO (1+1-RO (I) *KO (1) ) 

CONTINUE 

GO TO 7 

Dote ts), NCELL 

ETH1 (I) =2.4E-12*DN (I) * (2. 0*TEO (I-DTO (1I)) 

EK E1 (1) =2.091E-25*DN (I) * (UO (1) U0 (I+1) ) **2 

EB1(I) =B0 (I) *BO (I) /8.0* (RO (I1+1) *RO (1+1-RO (I) *RO (1) ) 

CONTINUE 

POWER CALCULATIONS. 

PABST=0. 

PBT=0.0 

PTHT=0,. 

PKET=0.0 

D0Se t= NCSLL 

PABS (I) =2.4E-19*DN (1) *HEAT (LI) 

PB(L) =1.9E-07* (ZB1 (I-EBO (I)) /DT 

PTH(1) =1.0E-07* (ETH1 (I-ETHO (I))/DT 

PKE (1) =1.0E-07* (EKE1(I-EKEO (I)) /DT 

PABST=PABST* PABS (L) 

PTHT=PTHT+PTH (1) 

PK ET=PKET+ PKE (1) 

PBT=PBT+t+PB(TI) 

CONTINUE 

ETHO,EKEQ,EBO ARE SET TO THE LAST VALUES, IE. ETH1,... 

DO4 I=1,NCELL 

ETHO (I) =ETH1 (1) 

EKEOQ (1) =EKE1 (1) 

EBO(1) =EB1(T) 

CONTINUE 

RETURN 

END 

; TPDATA 

sDAL#B RO (101) , UO (101) »TEO (100) ,DTC (190) ,NO(1C0) ,B0 (100) 

REAL VS (102) , TES (102) , TIS (102) ,NS (102) ,RS(102), BS (102) _ ,TPS (192) 
1,120 (100) ,TRF (100) ,TP2 (100) ,E3ZZ (100) , W341 (100) ,TRARD (100), 
2W3mM(100) 

EBON. PO, BZ 8 ,P)T]TE/COMHYO/U0, RO /COMTHO/TEO,DTO,NO,BO 


3 ARD,W3M 
1/PARAL/TPO TPF,TP2,£322, W3N1,TP ’ ; 
SUBROUTINE PREPARES NUMERICAL DATA TO BE STORED ON TAPE FOR 


——_> id 
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PLOTTING, 

DO1 I=1,P 

TES (1) =SNGL (TEOQ(I)) 

TIS (1) =SNGL (TEO(I-DTO (I) ) 

NS (I) =SNGL(NO (I) ) 

TPS (L) =TPF (1) 

CONTINUE 

N=P+1 

DOg2, c=1;N 

RS (1) =SNGL(RO(I)) 

BS (I) =SNGL (BO (I) ) 

VS (1) =SNGL (UO (I)) 

CONTINUE 

WRITECS} OT, IF, PO,R,BZ,P RS AQIS ARLES ISIN HINES Aas 
RETURN 

END 

SUBROUTINE FLUX(P,BO,80) 

REAL*8 BO (100) ,RO (101) 

INTEGER P 

TFLUX=0.0 

DO 1.-1=1,P 

TELUX=TFLOX +¢ BO (1) *3.14159* ((RO (1+1) **2) - (RO (1) **2)) 
CONTINUE 

WRITE (6,2) TFLUX 

FORMAT(* ',*TOTAL MAGNETIC FLUX IS',E14.6///) 
RETURN 

END 


This routine plots data from Magnetic tape of the M.H.D. routine. 


HiqQua 


REAL RM(102) , TE (102) ,TI (102) ,N(102),¥(102),B(102) ,TP (102) 
INTEGER P 

COMMON R,BZ,P0,T,TF,2,RM/PLOTTN/TE,TI,N,V,B ,TP 
TF IS TIME FOR FINAL PLOT. 

TINC IS INCREMENT IN TIME BETWEEN PLOTS 

TPL LS SET TO TIME FOR INITIAL PLOT. 

CALL PLOTS 

CAULEPEOT(2.0,5.0,-3) 

TF=40.0E-09 

TINC=5.0E-09 

TPL=20.CE-09 

READ (2, END=2)T,TF,P0,R,B2Z,2,RM,TE,TI,N,V,B ,TP 
IF(T.GE.TPL) GO TO 3 

GO TO 1 

TPL=TPLtTINC 

DO 4 'T=1,P 

TP (I) =Te (1) +T# (I) 

CALL HYDPLT 

IF(T.GE.TF) GO TO 2 

Go «rom 

CALL PLOT(0.0,0.0,999) 
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STOP 
END 

SUBKOUTINE HYDPLT 

SUBROUTINE CSALE FINDS SCALING PARAMETERS FOR COMPUTER PLOTTING 
AND RESULTS APPEAR IN LOCATION (P+1) AND (+2). 

AN ORDINATE OF 7 INCHES IN A PLOT IS ASSUMED, 

REAL TE (192) ,TI (102) ,N (102) ,B(102) ,V(102) ,RM (102) ,TP(102) 
INTEGER P,NLB(25) /68,69, 85,98, 73,99,104,64,77, 103, 113,112, 46,113 


1,120,47,47,67,84,46,96,115,47,47,93/,VLB(Z9) /70,83,86,102, 64,101 
2,69,83,86,67,73,99,104,64,77;61, 113,112, 46,118, 47, 47,67,84,97,98 
3709,01,93/7,VMAX 


COMMON R,B2,P9,T,TF,P,RM/PLOTTN/TE,TI,N,V,B 7mLe 

MXT,NOMAX,VMAX REPRESENT MAX. VALUES IN TEMP, DENSITY,AND VELOCITY 
TO A POWER OF 10 

ALL SCALING COEF. ARE CONTAINED IN LOCATIONS (P+1) AND (P#2). 
MAGNETIC FIELD SCALING- B/BZ K DENSITY (TIME=0)/HAX. DENSITY. 
PLOT WINDOW, X-AXIS AND LABEL PARAMETERS. 

CALL HYDLBL 

PLOT OF ELECTRON, ION, AND PARALLEL TEMPERATURE VS. RADIUS. 

CALL CSALE(TP,MXT, P,Q) 

TE (P+1) =TeP(P+1) 

TE (P+2) =TP (P+2) 

TI (P+1) =TE(P+1) 

TI (P+2) =TE(P +2) 

CALL AXIS(0.0,0.0,"TEMPERATURE (EV)*,16,7.0,90.0,TP(P+1) ,TPe (P+2) 
1,20.0) 

CALL DLINE(RM,TE,P,2,C.050) 

CALL DLINE(RM,TP,2,1,0.050) 


PLOT OF PLASMA DENSITY VS. RADIUS. 

CALL CSALE(N, NOMAX,P,1) 

CALL AXIS(7.5,0.0,' ',-1,7-0,90.0,N(2+1) ,N(P+2) ,26.0) 
NLB(15) =102+ NOMAX 

CALL LSYMB(7.9,2+5,0.12, NLB (1) ,90.0,-1) 

DO1 T=2,25 

CALL LSYMB(999.0,999.0,0.12,NLB(I) ,90.0,-1) 

CALL FLINE(RM,N,P,1,C,0) 


PLOT OF MAGNETIC FIELD VS. RADIUS. 
B(P+1) =N (P+1) 

B(P+2) =N(P+2) 

DO2 I=1,P 

B (I) =7.0E+16*PO*B (I) /(BZ*10.0**NOMAX) 0.005 
CALL DASHLN(RM,B,P,1) 


PLOT OF FLOW VELOCITY VS. RADIUS. 

WINDOW PLUS LABELLING, 

CALL PLOT(-0.5,10.0,-3) 

CALL HYDLBL 

COMPLETION OF GRAPH PLUS LABELLING. 

CALL SCALE(YV 72072 4,1) 

CALL AX1S(0.0,0.0,'FLOW VELOCITY (CM/SEC)',+22, 


17.0,90.0,V(P+1),V (P+2) ) 
CALL AXIS (76 5,0. 0, ‘FLOW VELOCITY (CM/SEC)',-22, 
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17.0,90.0,V(2+1),V (P+2)) 
CALL PLINE(RM,V,P,1,0,0) 


NOVE ORIGIN TO PREPARE FOR NEXT PLOT 
CALL PLOT (12.0,-11.0,-3) 

RETURN 

END 

SUBROUTINE HYDLBL 

REAL RM (102) 

INTEGER P 

COMMON R,BZ,PO0,T,TF,P,RM 


SUBROUTINS HYDLBL IS A LABELLING ROUTINE. GIVES A WINDOW (8 X 8.5), 


THE X-AXIS AND THE FOLLOWING PHYSICAL CONSTANTS, TIME(NSEC), 
FILLING PRESSURE (TORR), AND THE MAGNETIC FIELD (KG). 

R= MAX. VALUE FOR R (CM) 

T= TIME (SEC) 

DRAWING OF A WINDOW. 

CALL PLOT(0.0,0.0,3) 

GALL PLOT (0. 0,8.0, 2) 

CALL PLOT(8.5,8.0,2) 

GALE PLOY (8. 5;0.0, 2) 

CALL PLOT(0.0,0.0,2) 

DEFN OF A NEW ORIGIN. 

CALL PLOT(0.5,0.5,-3) 

SCALING COEF FOR RADIUS. 

RM (P+1) =0. 

RM (2+2) =R /7.0 
CALL AXIS(C.C,0.0,*RADIUS (CM) *,-11,7.5,0.0,RM (P+1), RM (P+2) ,14.29) 
CAUL PAXIS (0 07760 po fol lee Oe dp) VETO, 120,28. 57) 
PRINTING OF PLASMA PARAMETERS. 

BZKG=BZ¥*1.0E-03 

TNSEC=T*1.,0E+9 

CALL SYMBOI,(0.25,7.2,0.08,'TIME IS NSEC.*,0.0,19) 
CALL NUMBER (C.95,7.2,0.08,TNSEC,0.0,+2) 

CALL SYMBOL(2.1,7.2,0.08,'FILLING PRESSURE IS TORR. *,0.0,30) 
CALL NUMBER (3. 687,7.2,0. 08, 20,0.0,+1) 

CALL SYMBOL(4.85,7-2,0.08,"EXTERNAL FIELD IS KGAUSS.',0.0, 
130) 

CALL NUMBER (6.25,7.2,0.08,BZKG,0.0,1) 

RETURN 

END 

SUBROUTINE CSALE(V,MX, 2, SREQ) 

REAL V (102) 

INTEGER P,SREQ 

SUBROUTINE SCALES VARIABLES FOR PLOTTING. 

VALUES TO BE PLOTTED ARE NORMALIZED TO A POWER OF 10 (MX). 
SL=SCALE LENGTH/INCH 

SREQ= 0 (NO SCALING) =1 IF SCALING IS REQUIRED. 

SCALING PERFORMED ON BASIS OF 7 - INCH ORDINATE. 

X=0. 

p01 I=1,e 

IF(X.LT.V(I)) X=V (1) 

CONTINUE 

DIFF=ALOG10 (X) 

MX=IFIX (DIFF) 
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DIFF=DIFF-MX 
IF (DIFF.LT.0.813) G0 LORS 
MX=MX+41 
DIFF=DIFF-1 
a) DP USHEOQ. NEST) GOUTOTA 
SF=10.0**MX 
DO2 I=1,P 
V (I) =V (1) /SP 
CONTINUE 
LE(DIFF.GT.02512)6O.T0 5 
EP(DE'F.GT.0.2040)GO0 To 6 
SL=0.25 
lO) Yeo) 7/ 
5 SL=1.0 
GO TO 7 
SL=0..5 
IP (SKEQ.NE.1) SL=SL¥*10.0**MX 
V(P+1) =0.0 
V (P+2) =SL 
RETURN 
END 
SUBROUTINE DLINE(X,Y,N,ICHAR,HT) 
REAL X (102) ,¥(102) ,X1(102),¥1(102) 
SANDYS ARE VECTORS = TOMBE LOLRE Ds. SCALING MUST 
N=NO. OF PTS. TO BE PLOTTED. ICHAR=SYMBOL CODE. 
BOe1L 14 1,.N 
X17 (1) =(X (1-X (N+1)) /X (N+2) 
1 ¥1(I) =(Y(I-Y¥ (N#1)) /¥ (N+2) 
CALL SYMBOL (X7¢1)',¥101) (HT, ICHAR, 0.03 -7) 
DO 2 I=2,N 
2 CALL SYMBOL (XJ (1) pY1 (1), HT; 1CHAR 0.0, -2) 
RETURN 
END 
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BE ANNOTATED 
HT=CHARACTER HEIGHT, 
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This routine calculates a grid of (K**2*D]. 
Method of calculation depends on region in complex plane, The three 
methods are: 1. Integration as calculated by J. McMullin; 


DER 


3 
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. Power series expansion as given on page 178 of Stix. 


COMPLEX*8 F,CMPLX 
REAL*8 A,B,WR,WI,K,DSQRT 
REAL FPL? (101,101) ,EPLT1 (5050) ,FPLT2 (5151) 
EQUIVALENCE(FPLT (1,1) ,FPLT1) , (FPLT (1,51) ,FPLT2) 
K=2.0D0 
WMAX=10.0 
BASE=WMAX/100.0 
WRITE(7,8 
8 SP eTCETERA WR/WPE WI/WPE ALPHA 
1 FUNCTION ABS 1/ABS"//) 


pomeri=1,104 


Asymptotic expansion as given on page 180 of Stix, or 


BETA 


os Sb ey tasy aera 
warn Dea Ub wkusd 


a ats oe 
Seder ax 
bea oro ; 

a oe 
Bey eer a 
nas) Ce foe { i” i rit 


749 AK re Ta | sad 
Bren a" 


ie “@oe7e 
2£2t,, ta oe, 
4% 9 
Prey ve mel 
At. u. RaSea 


ne 


= ND 


pOr2 12=1;,101 

WR=BASE* (I1-1) 

WI=BASE* (12-1) 

A=WR/ (1.4142 13562*kK) 

B=WL/ (1.414213562*K) 

XRAD=DSQRT (A**2+B**2) 

BS=B 

IF(XRAD.LT.2.8) GG TO 3 

IF(BSaLT.9,2)) GO TO. 4 

CALL INTG(A,B, F,K) 

GOTO 5 

CALL EXPAN (A,B,F,K) 

GOM TO 5 

CALL ASYMP(A,B,F,K) 

CONTINUE 

FPLT (12,11) =CABS (F) 

XINF=1.0/FPLT (12,11) 

TF ((( (12/2) *2) .NE.I2) .AND. (( (11/2) *2) .NE.11)) 
TWRITE (7,7) WR,WI,A,B,F,FPLT (12,11) ,XINF 

PORMADAt ie 2 DAGser 1X5010.3, 3X), 2 (E10. 3, 1k R10. asak)) 
CONTINUE ; 

CONTINUE 

WRITE(4) FPLT1 

WRITE(4) PPLT2 

STOP 

END 

SUBROUTINE INTG(A,B,F,K) 

COMPLEX*8 EF, CMPLX 

RE AL*8 A,B, X (10) -W(10) , INT1T,IN&2,F1,F&2,F 1M, 720,48, ABZ, A2B2 
a, KsEuer Lp DEXP, DCOS,DSIN 

DATA X/C0. 245 3407983009D0, 0. 7374737285454 DC, 1. 2340 762153953D0 
1,1.7385377121166D0 ,2.2549740020893D0, 2. 78880 6058428106 
2,3.347854567 38 3200, 3. 944764040 1156D0,4 .6036824495507 DO 
3,5.2874808900112D0/, 
UW/4. 6224 36696006D-1, 2. 866755053628D-1,1.090172060200D-1 
5, 26481052038 746D-2, 3. 243773 342238D-3, 2.283386360163D-4 
6,7. 802556478532D-6,1.086069370769D-7,4.399340992273D-10 
7726229393645 534D-13/ 

INT1=0.0D9 

INT2=0.0D0 

Doel i=1, 10 

F1=X (1) /( (X (I-A) **2+B¥*2) *0. 5641895835 

F2=F1*X (1) 

XM=-X (I) 

FIM=XM/ ((XM-A) **24+B¥*2) *0.5641895835 

F2M=FIM*XM 

INTI=INT1+W (1) *(F1+F1M) 

INT2=INT2+wW (1) * (F2+F 2M) 

CONTINUE 

AB2=2.0*A*B 

AZ B2=A* *Z-B*¥ *2 

IF Cie ony 

AZB2.LT.-150.0D0) A2B2=- : 
ae nae e 3 544907702*DEXP (-A2B2) * (BYDCOS (AB2-A*DSIN (ABZ) ) #1NT2 
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PI=3.544907702*DEXP (-A2B2) * (A*DCOS AB2) + B*¥D B2)-B* 
F=CMPLX(SNGL (FR) ,SNGL (FL) ) poe eye eee 
RETURN 

END 

SUBROUTINE EXPAN (A,B,F,K) 

COMPLEX*16 2,DCMPLX,S,S1,F0,J,CDEXP 

COMPLEX*8 PF 

REAL*8 K,A,B,DIV,CDABS 

J=DCMPLX(0.0D0,1.0D0) 

Z=DCMNPLX (A, (-3B)) 

DIV=1.0D0 

S=DCMPLX (0.0D0,0.0D0) 

IT=0 

etal 

Semen RUNES eee Tietze U00 eee oy 
S=St 

DIV=DIV* (2.0D0*1+1.0D0) 

IF (SNGL (CDABS (S1)) «GT. SNGL(1.C0D-3*CDABS(S))) GO TO 1 
PO=1.772453851*CDEXP (- (Z¥**2) ) +2. 0D0*J*S 
F=K*¥*24+1,0D0 +3 *Z*FO 

RETURN 

END 

SUBROUTINE ASYMP (A,B,F,K) 

COMPLEX*16 Z2,DCMPLX,J,S,CDEXP, FO 

COMPLEX*8 F 

REAL*8 A,B,K 

J=DCMPLX(0.0D0,1.0D0) 

Z=DCMPLX (A, (-B)) 

S=0.5D0/Z+0. Z5D0/(Z2**3) +0.375D0/ (2**5) +0.9375D0/ (2**7) 
1+3.28125D0/(2**9) 

FO=1. 7724538 51*CDEXP (> (Z**2) )+2.0D0¥*J*S 
F=K**2+1.,0D0 +0 *Z*FO 

RETURN 

END 


This routine makes a contour map of data from previous routine, 
CCONTR is a subroutine in *CSLIB at computing center, U. of A. 


REAL A(191,101),CVAL(10) ,AR1 (5050) ,AR2(5151) ,MINVAL 
INTEGER IOP(4) 
DATA I0P/1,1,0,0/ 
EQUIVALENCE (A(1,1) ,AR1), (A(1,51) ,AB2) 
READ(4) AR1 
READ(4) AR2 
MINVAL=1.0 
DO 1 I=1,10 
CVAL (I) =I*MINVAL 
1 CONTINUE 
XSIZE=10.0 
YSIZE=10.0 
M=101 
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CALL PLOTS 

Cc CALL VRTWIN(0.0,13.78,0.0,10.5) 
CALLS PLOT (S.0;5.0, -3) 
CALL CCONTR(XSIZE, YSIZE,A,M,N, CVAL, NC, IOP) 
CALL PLOT(0.0,0.0,999) 
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This routine makes a 3-D plot of data from previous routine. 
CPLT3D is a subroutine in *CSLIB at computing center, U. of A. 


REAL A(191,1C1),AR1 (5050) ,AR2(5151) ,MAXVAL 
EQUIVALENCE(A(1,1) ,AR1), (A(1,51) , AR2) 
READ(4) AR1 
READ(4) AR2 
MAXVAL=1.0 
DOwtal=1, 101 
DO 2 J=1,101 
A(J,1I) =1.0/A (J,1I) 
IF(A(J,1).GT.MAXVAL) A(J,1) =MAXVAL 

2 CONTINUE 

1 CONTINUE ; 
M=101 
N=101 
CALL PLOTS 

c CALL VRIWIN (050,13.78,0.0,10.5) 
SCALine PLOTS. 0p. 0, -3) 

K=3 
DIST=5.0 
YAW=330.0 
PITCH=30.0 
SIZE=10.0 
KODE=0 
MGN=2 
SCALE=1.0 Fe Jam 
CALL CPLT3D(A,N,M,K,DIST,YAW, PITCH, SIZE, KODE, NGE,SCALE) 
GALDMELOT (0, 0; 0. 0, 999) 
STOP 
END 


-<_=- ww" ee ew ew eo oe — ——- osname Ganesan adhesin noe aes == aS === ~~ 
“<< - - 
=_- 


This routine calculates total energy in plasma from tape data. 


45 (102), BS (102) 

REAL VS (102) , TES (102) ,TIS (102) , NS (102) , RES (1 

1, TPS (102) ,RD (102) , EB (102) ,EK (102) ,ET (102) ENT (102) 
INTEGER P,PM1- 
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THIS ROUTINE READS MAG, TAPE AND CALCULATES TOTAL PLASMA ENERGY, 


WRITE (5,1) 
FORMAT(* '," ENTER IX=1 IF TAPE CONTAINS INFO. FROM NONLINEAR HEATI 


ING, IX=2 1F INV. BREM, ONLY. ENTER TAGE RECORD NO, '/ 
Zien TX. CREECH) 


READ (5,2) IX,NREC 
FORMAT (215) 
IZ=NREC-1 
TR iZe to. 0) CO<TOES 
DO 12 1=1,12 
READ (3) 
CONTINUE 
CONTINUE 
GOeTO (4,5), 7x 
READ(3)T,TF,P0,R,B2,P,RMS, TES, TLS,NS,VS,BS,TPS 
On 1040 T=), 
TES (1) =TES (1) #0. 33333*TPS (L) 
GO TO 6 
READ(3)T,TF,PO,R,B2Z,P, RUS, TES, TIS,NS,VS,BS 
CONTINUE 
EBT=0.0 
EKT=0.0 
ETT=0.0 
DO 8 I=1,P 
VS (I) =0.5*(VS(1) #VS(I+1) ) 
EB (1) =0.125* (BS (1) **2) * (RNS (I+ 1) **2-RBMS (1) **2) 
EK (I) =2.6ZE-24¥*NS (L) * (VS (1) **2) * (BMS (I+1) **2-RMS (1) **2) 
ET (1) =7. 54E-12¥*NS (1) * (TES (1) +TIS (1) ) * (RMS (I+1) **2-RMS (I) **2) 
EBT=EBT+EB (1) 
EKT=EKT+EK (1) 
EVT=ETT+ET (I) 
CONTINUE 
ENG=LBT+EKT+ETT 
WRITE(6,20)T 
FORMAT('0*,*TIME IS',E12.3,' SECS."') 
PE ENG, ETT,EBT,EKT 
EMI iat ee acieke ENERGY = ',E12.5/" ',' THERMAL ENERGY = ',E12.5/ 
' * / 'MAGNETIC ENERGY = ',E12.5/" *,*KINETIC ENERGY = ',E12,5) 
pO 1¢-I=1,e 
ENT (1) =EB (I) +EK (1) +ET (I) 
WRITE (6,11) ENT(I) ,ET (I) ,EK (I) ,EB(Z) , RMS (I) 
FORMAT(* *,6E15.3) 
CONTINUE 
STOP 
END 
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